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Abstract. The problem of slowing down light by orders of magnitude has been 
extensively discussed in the literature. Such a possibility can be useful in a variety 
of optical and microwave applications. Many qualitatively different approaches 
have been explored. Here we discuss how this goal can be achieved in linear 
dispersive media, such as photonic crystals. The existence of slowly propagating 
electromagnetic waves in photonic crystals is quite obvious and well known. The 
main problem, though, has been how to convert the input radiation into the 
slow mode without losing a significant portion of the incident light energy to 
absorption, reflection, etc. We show that the so-called frozen mode regime offers 
a unique solution to the above problem. Under the frozen mode regime, the 
incident light enters the photonic crystal with little reflection and, subsequently, 
is completely converted into the frozen mode with huge amplitude and almost 
zero group velocity. The linearity of the above effect allows the slowing of light 
regardless of its intensity. An additional advantage of photonic crystals over other 
methods of slowing down light is that photonic crystals can preserve both time 
and space coherence of the input electromagnetic wave. 
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1. Introduction 

1.1. What is slow light? 

It is common knowledge that, in vacuum, light propagates with constant velocity 
c w 3 x 10 8 m/sec. In optically transparent nondispersive media, the speed of light 
propagation is different 

v = uj/k = c/n, (1) 

where k is the wave number, ui is the respective frequency, and n is the refractive 
index of the medium. At optical frequencies, the refractive index n of transparent 
materials usually does not exceed several units, and the speed of light propagation is 
of the same order of magnitude as in vacuum. 

The situation can change dramatically in strongly dispersive media. Although 
the phase velocity of light is still determined by the same expression Q, the speed of 
electromagnetic pulse propagation is different from v and is determined by the group 
velocity PEHS] 
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which is one of the most important electromagnetic characteristics of the medium. 
With certain reservations, the group velocity u coincides with the electromagnetic 
energy velocity and is usually referred to simply as the propagation speed of light in 
the medium. Hereinafter, the speed of light propagation means the group velocity J2J), 
rather then the phase velocity QJ. 



Slow light in photonic crystals 



2 



Strong dispersion means that the group velocity u strongly depends on the 
frequency and can be substantially different from c. In the slow light case, which is the 
subject of our interest, the electromagnetic pulse propagates through the dispersive 
medium at the speed u <C c, regardless of the respective value of the phase velocity Q . 
In some cases, u can even become vanishingly small implying that the propagating 
electromagnetic mode at the respective frequency does not transfer energy. In another 
extreme case, the group velocity u can exceed c (the so-called case of superluminal 
pulse propagation) , without contradicting the causality principle P El 13 ■ In yet 
another case of a left-handed medium, the group velocity u can have the opposite sign 
to that of the phase velocity v |SJ. But again, in this paper we will focus exclusively 
on the slow light and related phenomena. 

Slow and ultraslow light have numerous and diverse practical applications. The 
related phenomena include dramatic enhancement of various light-matter interactions 
such as nonlinear effects (higher harmonic generation, wave mixing, etc.), magnetic 
Faraday rotation, as well as many other important electromagnetic properties of the 
optical media. Such an enhancement can facilitate design of controllable optical 
delay lines, phase shifters, miniature and efficient optical amplifiers and lasers, etc. 
In addition, ultraslow light might allow nonlinear interactions down to a single 
photon level, which could significantly benefit the design of ultrasensitive optical 
switches, quantum all-optical data storage and data processing devices. Ultraslow 
light can also be used in quantum communication and design of novel acousto- 
optical devices. This list can be continued. For more detailed information on 
the prospective practical applications of slow light phenomena see, for example, 
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and references therein. 

1.2. Temporal dispersion versus spatial dispersion 

In recent years, several different approaches have been pursued in order to slow down 
or even completely stop light. These approaches can be grouped into two major 
categories: 

- those where the low group velocity results from strong temporal dispersion dn/duj 
of optical media; 

- those where the low speed of pulse propagation is a result of coherent interference 
in spatially periodic heterogeneous media, such as photonic crystals. 

Let us start with a brief discussion of slow light phenomena in media with strong 
temporal dispersion. 

1.2.1. Slow light in media with strong temporal dispersion Assuming that the 
refractive index n in J5J is of the order of unity, which is usually the case at optical 
frequencies, one can state that a very low group velocity can only occur if n varies 
strongly as a function of u 



Strong frequency dependence of the refractive index n can be a result of excitation 
of electronic or some other intrinsic resonances of the medium, which are normally 
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accompanied by strong absorption of light. Recently, though, several techniques have 
been developed that allow to significantly suppress the absorption of light at the 
frequency where the derivative dn/duj peaks. 

One of the most successful ways to suppress light absorption is based on the effect 
of electromagnetically induced transparency (EIT) 125]. In such a case, the incident 
light interacts with atomic spin excitations forming combined excitations of photons 
and spins, called dark-state polaritons. These polaritons propagate slowly through 
the medium in the form of a sharply compressed pulse, the energy of which is much 
smaller than that of the incident light pulse. Most of the incident light energy is 
expended to create the coherent state of the atomic spins, the rest is carried away 
by the control electromagnetic field. The pulse delay inside the medium is limited 
by the bandwidth of the transparency window, which decreases with propagation 
distance. At higher propagation distances the medium becomes increasingly opaque 
at frequencies other than the line center, further reducing the available transparency 
window |12II13) . Specific physical mechanisms of such transformations are very diverse. 
The detailed description of EIT and related phenomena can be found in the extensive 
literature on the subject (see, for example, EH EU El E3 E], and references 
therein). The techniques based on EIT have already produced some amazing results, 
such as reduction of the speed of pulse propagation by 7 - 8 orders of magnitude 
compared to the speed of light in vacuum, while providing a huge and controllable 
pulse delay. 

Another method to create a transparency window in otherwise opaque substance 
was used in QUI • This method involves the creation of a spectral hole by the 
periodic modulation of the ground state population at the beat frequency between 
the pump and the probe fields applied to the material sample. It can produce slow 
light in a solid-state material at room temperature. The spectral hole created by this 
technique can be extremely narrow (36 Hz in the experiment |151 116p . and leads to 
a rapid spectral variation of refractive index. It allowed to reduce the light group 
velocity in a ruby crystal down to 57 m/s. 

Physically, the above approaches to slowing down the light can be viewed as 
a reversible transformation of the input nearly monochromatic light into some kind 
of coherent atomic excitations (e.g., dark-state polaritons) with very low relaxation 
rate and low group velocity. In other words, ultraslow pulse propagating through 
such a medium is, in fact, an intrinsic coherent excitation triggered by the input 
light, rather than a light pulse per se. This process always involves some kind of a 
delicate resonant light-matter interaction with extremely small bandwidth. Indeed, 
the relation JSJ yields the following limitation on the slow pulse bandwidth 



where the assumption is made that the refractive index n within the transparency 
window is of the order of unity. The condition (@J can also be viewed as a constraint 
on the minimal propagation speed of a light pulse with a given bandwidth Acj. On 
the positive side, the approach based on EIT or its modifications does produce an 
exceptionally low speed of pulse propagation, which can have some very important 
practical implications. 

In the rest of the paper we focus exclusively on those techniques which do 
not involve any intrinsic resonant excitations of the medium and, therefore, do not 
essentially rely on strong temporal dispersion. Instead, we will focus on spatially 



Slow light in photonic crystals 



4 



periodic dielectric arrays, in which low group velocity results solely from spatial 
inhomogeneity of the optical medium. 

1.2.2. Slow light in spatially periodic arrays Well-known examples of optical periodic 
dielectric structures include photonic crystals periodic arrays of coupled optical 
resonators [HO d 1201 El 1221 031 121 , and line-defect waveguides [23]. Generally, a 
periodic heterogeneous medium can be assigned a meaningful refractive index n only 
if the structural period L is much smaller than the light wavelength A 

L«A. (5) 

On the other hand, a significant spatial dispersion associated with heterogeneity of 
the medium can occur only when L and A are comparable in value 

L ~ A. (6) 

In particular, the relation 10 defines a necessary condition under which heterogeneity 
of the medium can lead to low speed of electromagnetic pulse propagation. Hence, in 
the cases where low speed of pulse propagation is a result of strong spatial dispersion, 
one cannot assign a meaningful refractive index to the composite medium, and the 
expression (J2Jl for the group velocity of light does not apply. 

At optical frequencies, the speed of pulse propagation in periodic dielectric arrays 
can be reduced by two or three orders of magnitude. This is not a fundamental 
restriction, but rather a technological limitation related to the difficulty of building 
flawless periodic arrays at nanoscales. On the positive side, the dielectric components 
of the periodic array are not required to display strong temporal dispersion and, 
hence, absorption of light is not an essential and unavoidable problem in this case. 
In addition, the photonic crystal based approach is much more versatile in terms of 
the input light intensity. It allows the same photonic device to operate both at high 
and low intensity of the input light. By contrast, utilizing strong temporal dispersion 
always involves significant nonlincarity and usually is limited to a certain amplitude 
of the input light. 

There is a natural bandwidth limitation on the slowed pulse in periodic dielectric 
arrays, which is similar to the case of slow light in time-dispersive media. Indeed, let 
Auj be the frequency bandwidth of a pulse and Ak - the respective range of the Bloch 
wave number. The average group velocity (u) of the pulse is defined as 

*■>"£■ < 7 » 

Let us make the following natural assumptions. 

(i) The pulse propagating inside the periodic medium is composed of the Bloch 
eigenmodes belonging to the same spectral branch of the dispersion relation lo (k) . 
This assumption implies that Afc cannot exceed the size 2n/L of the Brillouin 
zone 

Ak < 27T/L, (8) 
where L is the unit cell length of the periodic array. 

(ii) The refractive index of the constitutive components of the periodic array is of the 
order of unity and, therefore, 

L ~ A = 2ttc/w, (9) 
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where Ao is the light wavelength in vacuum. 

The relations Ij7l9|) yield the following limitation on the minimal speed of pulse 
propagation for a pulse with a given bandwidth Au 



The restriction (|10|l is similar to that defined by the inequality (0} and related to the 
case of slow light in a uniform medium with strong temporal dispersion. In either 
case, a higher refractive index would lower the minimal speed (u) of pulse propagation 
for a given pulse bandwidth Alo. 

Any attempt to circumvent the restriction (|10() would involve some kind of pulse 
compression techniques |17| . 

1.2.3. Examples of periodic arrays supporting slow light 

Coupled resonator optical waveguide. During the last several years, a tremendous 
progress has been made in theory and applications of periodic arrays of coupled optical 
resonators. Generally, if the coupling between adjacent resonators in a periodic chain 
is weak, the group velocity of Bloch excitations supported by such a periodic array is 
low. This is true regardless of the nature of individual resonators. The above simple 
idea forms the basis for one of the most popular approaches to slowing down the light. 
An extensive discussion on the subject and numerous examples and references can be 

found in [13 d 113 EH E2 E3 El E2 

A qualitatively similar situation occurs in line-defect waveguides in a photonic 
crystal slab, where a periodic array of structural defects plays the role of weakly 
coupled optical resonators. Following [23, consider a dielectric slab with a two- 
dimensional periodic array of holes in it. A row of missing holes in this periodic 
array forms a line defect, which supports a waveguiding mode with two types of cutoff 
within the photonic band gap. These characteristics can be tuned by controlling the 
defect width. Theoretical calculations supported by interference measurements show 
that the single waveguiding mode of the line-defect waveguide displays extraordinarily 
large group dispersion. In some instances, the corresponding traveling speed is 2 orders 
of magnitude slower than that in air. According to |25| . one of the major limiting factor 
here is structural imperfection of the array. 

Slow light phenomena in periodic arrays of weakly coupled resonators have been 
the subject of a great number of recent publications, including some excellent review 
articles cited above. For this reason, further in this paper we will not discuss this 
subject any more. 

Photonic crystals. Photonic crystals are spatially periodic structures composed of 
usually two different transparent dielectric materials jlffij . Similarly to periodic arrays 
of coupled resonators, in photonic crystals, a low group velocity of light can result 
from multiple scattering of individual photons by periodic spatial inhomogeneities, 
rather than from temporal dispersion of the substance [22 H3 E3 EH1 El CM EI] The 
lowest group velocity achievable in photonic crystals for a given pulse bandwidth can 
be close to that defined by the fundamental restriction (|1C)|> . For example, if we want 
a pulse to propagate undistorted at speed as low as 10 -3 c, its bandwidth Aw should 
be less than 10 _3 w, which at optical frequencies is of the order of 10 GHz. In this 
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respect, the situation in photonic crystals is as good as it can possibly be in any other 
linear passive media with limited refractive index. 

Unlike the case of optical waveguides and linear arrays of coupled resonators, 
in photonic crystals we have bulk electromagnetic waves capable of propagating in 
any direction through the periodic heterogeneous structure. This results in much 
greater density of modes, compared to that of the above-mentioned arrays of coupled 
resonators. In addition, electromagnetic waves in photonic crystals can remain 
coherent in all three dimensions, which is also essential for a variety of practical 
applications. 

A major problem with slow light in photonic crystals is the efficiency of conversion 
of the incident light into the slow mode inside the heterogeneous medium. We shall 
see in the next section that in most cases an incident electromagnetic wave with the 
frequency of one of the slow modes is simply reflected back to space, without creating 
the slow mode inside the photonic crystal. How to overcome this fundamental problem 
and, thereby, how to transform a significant fraction of the incident light energy into 
a slow mode with drastically enhanced amplitude, is one of the primary subjects of 
this paper. 

The paper is organized as follows. In Section 2 we describe, in general terms, 
what kind of slow modes can exist in photonic crystals and under what circumstances 
some of these modes can be effectively excited by incident light. We show, that there 
is a unique situation, which we call the frozen mode regime, in which the incident light 
can enter the photonic crystal with little reflection and be completely converted into 
a slow mode with nearly zero group velocity and drastically enhanced amplitude. 

Section 3 gives an overall picture of the frozen mode regime in periodic layered 
media, without going into the detailed analysis based on the Maxwell equations. All 
the statements made in this section are later proven in Sections 5 through 11. 

In section 4 we define the physical conditions under which a periodic layered 
array can support the frozen mode regime. These conditions boil down to whether 
or not the electromagnetic dispersion relation of the periodic array can develop a 
stationary inflection point (|15fl . This requirement imposes quite severe restrictions on 
composition and geometry of the periodic layered medium. We show, in particular, 
that in the case of light propagating normally to the layers, the frozen mode regime can 
only occur if some of the layers are magnetic with significant nonreciprocal Faraday 
rotation. In the case of oblique light propagation, the presence of magnetic layers 
is not required, which makes it possible to realize the frozen mode regime at any 
frequency range, including optical and UV. A trade-off though is that at least some of 
the layers of a non-magnetic stack must display significant dielectric anisotropy with 
tilted orientation of the anisotropy axis. 

Section 5 is devoted to electrodynamics of periodic layered media. Particular 
attention is given to the cases where some of the layers display dielectric and/or 
magnetic anisotropy, because otherwise, the electromagnetic dispersion relation u> (k) 
of the periodic array cannot develop a stationary inflection point (|15fl and, therefore, 
such a structure cannot support the frozen mode regime. 

Sections 6 through 12 constitute the analytical basis for the entire investigation. 
There we present a rigorous and systematic analysis of the scattering problem for a 
semi-infinite periodic array of anisotropic dielectric layers. The emphasis is on the 
vicinity of stationary points (|12fl of the electromagnetic dispersion relation, where 
the slow electromagnetic modes can be excited. The comparative analysis of all 
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possible stationary points shows that only a stationary inflection point (|15f) can provide 
necessary conditions for slowing down and freezing a significant fraction of incoming 
radiation. In all other cases, the incident wave is either reflected back to space, or 
gets converted into a fast propagating mode with low amplitude. The exact analytical 
results of these sections are supported by a number of numerical simulations. 

2. Stationary points of dispersion relations and slow modes 

In periodic heterogeneous media, such as photonic crystals, the velocity of light is 
defined as the wave group velocity 

u^duj/dk, (11) 

where k is the Bloch wave vector and uo — ui (k\ is the respective frequency. At some 
frequencies, the dispersion relation uj ( k J can develop stationary points 



du/dk = 0, (12) 

where the group velocity u vanishes. Zero group velocity usually implies that the 
respective Bloch eigenmode does not transfer electromagnetic energy. Indeed, with 
certain reservations, the energy flux S of a propagating Bloch mode is 

S = Wu, (13) 

where W is the electromagnetic energy density associated with this mode. If W is 
bounded, then the group velocity u and the energy flux S vanish simultaneously at the 
respective stationary point (|12|) of the dispersion relation. Such modes are referred to 
as slow modes, or slow light. Some examples of stationary points Ijl2(l are shown in 
Fig. ^ where each of the frequencies oj a , u>b, ojq is associated with a slow mode. 

The electromagnetic dispersion relation of any photonic crystal displays an infinite 
number of stationary points like those shown in Fig. ^ But, a common problem with 
almost all of them is that the respective slow modes cannot be excited in a semi- 
infinite photonic crystal by incident light. This explains why there have been only a 
few attempts to exploit the photonic crystals for slowing down the light. Let us take 
a closer look at this problem. 

Consider a scattering problem of a plane monochromatic wave normally incident 
on a lossless semi-infinite photonic slab with the electromagnetic dispersion relation 
shown in Fig. ^ The symbol k in Fig. denotes the normal component of the Bloch 
wave number k in the periodic structure, which in the case of a normal incidence is 
the only nonzero component of k. The symbols ^r, and ^t in Fig. [21 denote 
the incident, reflected, and transmitted waves, respectively. The transmittance r and 
reflectance p of the semi-infinite slab are defined as 

St S R 

T = -a", P= FT" = 1-T. (14) 

where Si, Sr and St are the normal energy fluxes of the respective waves. 

If the frequency oj is close to the band edge frequency uj g in Fig. ^ then the 
incident wave will be totally reflected back into space, as illustrated in Fig. |3| This 




Figure 1. An example of electromagnetic dispersion relation u)(k) with various 
stationary points: (i) extreme points a and b of the respective spectral branches, 
(ii) a photonic band edge g, (iii) a stationary inflection point 0. Each stationary 
point is associated with slow light. 



implies that the fraction of the incident wave energy converted into the slow mode 
corresponding to the point g in Fig. ^ vanishes as u — > ui g . 

In another case, where the incident wave frequency is close to either of the 
characteristic values to a or w;, in Fig. ^ the slab transmittance remains finite, as 
seen in Fig. [3J This implies that the incident wave will be partially transmitted into 
the semi-infinite photonic slab. The problem, though, is that none of the transmitted 
light will propagate inside the slab in the form of the slow mode corresponding to the 
respective stationary point a or b. For example, at frequency u> a , all the transmitted 
light corresponds to a fast propagating mode with positive and large group velocity 
and the wave number different from that corresponding to the point a in Fig. ^ A 
similar situation takes place at u> = the fraction of the transmitted light that is 
converted into the respective slow mode vanishes as oj — > u>b- 

Let us turn now to the stationary inflection point in Fig. where both the 
first and the second derivatives of the frequency oj with respect to k vanish, while the 
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Figure 2. Plane wave normally incident on a lossless semi-infinite photonic slab. 
The subscripts I, R, and T refer to the incident, reflected and transmitted waves, 
respectively. 



third derivative is finite 

dui d 2 u> d^oj 

at uj = uj and k = k : — 0; = 0; > 0. (15) 

In such plane wave with lo = loq incident from the left can be transmitted 

into the semi-infinite photonic crystal with little reflection, as demonstrated in Fig. 
121 But most remarkably, having entered the photonic slab, the light is completely 
converted into the slow mode with infinitesimal group velocity and drastically 
enhanced amplitude. Such a behavior is uniquely associated with stationary inflection 
point (|15|l of the dispersion relation and constitutes the frozen mode regime BTfllHHl BT] . 
In the frozen mode regime, the vanishingly small group velocity u in Eq. Ijl3(l is offset 
by the diverging value of the energy density W 

as to — * loq: u ~ \u> — luo\ 2 ^ 3 — > 0, W ~ \lu — luq\ 2 ^ 3 — > oo, (16) 

As a result, the energy flux (|13|) associated with the transmitted frozen mode remains 
finite and comparable with that of the incident wave even at the frozen mode frequency 
ujq corresponding to the point of the dispersion relation in Fig. ^ Such a 
spectacular behavior is uniquely attributed to a stationary inflection point H15|l of 
the electromagnetic dispersion relation. Of course, in reality, the electromagnetic 
energy density W of the frozen mode will be limited by such factors as absorption, 
nonlinear effects, imperfection of the periodic dielectric array, deviation of the incident 
radiation from a perfect plane monochromatic wave, finiteness of the photonic slab 
dimensions, etc. Still, with all these limitations in place, the frozen mode regime can 
be very attractive for a variety of practical applications. 

In the following sections we present a detailed analysis of the frozen mode regime 
associated with stationary inflection point Q15[l. In the rest of this section we briefly 
discuss the effect of photonic crystal boundaries on slow light phenomena. 
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Figure 3. Transmittance r of the semi-infinite photonic slab as a function of 
incident light frequency ui for the semi-infinite photonic slab with the dispersion 
relation presented in Fig. 1. The characteristic frequencies ui a , Wj, ujq, and w g are 
associated with the respective stationary points in Fig. 1. Within the photonic 
band gap at ui > ui g the incident light is totally reflected by the slab. 



2.1. Slow light in a finite photonic slab 

Up to this point we have considered light incident on the surface of a semi-infinite 
photonic crystal. Since real photonic crystals are always bounded, the question arises 
whether and how the photonic crystal boundaries affect the conditions of slow mode 
excitation and propagation. 

To start with, let us recall that in an unbounded (infinite) photonic crystal, the 
speed of light propagation is defined as its group velocity Ullfl. which determines the 
speed of pulse propagation in the medium. The spatial length I of a pulse inside the 
unbounded periodic medium is 

ii 

l~lo- (17) 

c 

where is the spatial length of the same pulse in vacuum. The quantity Iq is directly 
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Photonic Slab 




Figure 4. A pulse of length irj approaching a photonic slab of thickness D. The 
arrow shows the direction of pulse propagation. What happens after the pulse 
hits the slab boundary is shown in Fig. 5. 



Photonic Slab 



1 



D 



Figure 5. After hitting the slab, the pulse splits into the reflected and 
transmitted pulses. In a slow light situation, the transmitted pulse gets 
compressed in space. 



related to the pulse bandwidth Aw 

Au Xq 2t: c 



, , (18) 

IjJ Iq UJ Iq 



where 



A = — c 



is the light wavelength in vacuum. 

If instead of an infinite photonic crystal we have a bounded photonic slab of 
thickness D, as shown in Fig. 0]andEl the simple interpretation of the group velocity 
u as the speed of pulse propagation can still apply, provided that the pulse length I 
inside the photonic slab is much smaller than the slab itself 

(19) 

In other words, one can introduce the speed of pulse propagation inside the slab only 
if the entire pulse can fit inside the slab, as in the situation shown in Fig. \5\ In the 
slow light case, the group velocity u decreases sharply, and so does the pulse length I 
in (|17|l . Therefore, a slow pulse with a fixed bandwidth Alo is more likely to fit inside 
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the photonic slab than a fast pulse with the same bandwidth. The slower the pulse is, 
the better the condition l|19[) is satisfied. Taking into account the relations (|17f) and 
l|18fl . the condition l|19|) can also be recast as a lower limit on the pulse bandwidth 

2nu 

Aw>— , (20) 

implying that in order to fit inside the slab, the pulse bandwidth should not be too 
narrow. 

If a pulse satisfying the condition l|19|l or, equivalently, 1|20|) is incident on a finite 
photonic slab, the slab can be treated as a semi-infinite medium until the pulse actually 
hits the opposite boundary of the slab. Except for the next subsection, all the results 
discussed in this paper relate to the case (|19fl . where we can explicitly and literally 
talk about pulse propagation inside the medium and where the group velocity u in 
Ijllfl does have the meaning of the speed of pulse propagation. 

2.2. Resonance effects in a finite photonic slab 

A qualitatively different picture emerges if the pulse length I defined in Eq. I|17fl is 
comparable in magnitude or exceeds the slab thickness D. In such a case, the slab is 
too thin to accommodate the entire pulse and the electromagnetic field 'J'^ inside the 
slab becomes a superposition of forward and backward propagating waves undergoing 
multiple reflections from two opposite boundaries of the slab. This situation by no 
means can be interpreted as an individual pulse propagating through the periodic 
medium, because at any moment of time the electromagnetic field inside the slab 
cannot be viewed as a wave packet built around a single propagating mode. The term 
slow light does not literary apply here and, therefore, this case goes beyond the scope 
of this paper. Yet, it would be appropriate to discuss briefly what happens if the 
photonic slab becomes too thin to be treated as semi-infinite. 

Assume that the photonic slab is thin enough to satisfy the inequality 

Z>D, (21) 

which is opposite to (|19|1 . The condition H21(l establishes an upper limit on the incident 
pulse bandwidth 

27TU . . 

Aw « — . (22) 

Consider a plane monochromatic wave incident on a finite photonic slab in Fig. 
Since a monochromatic wave packet has I — > oo, the relations (|21|l and (|22|l are 
perfectly satisfied. If the photonic slab is lossless, its steady-state transmittance and 
reflectance are defined by the following expressions 

Sp S T Sr 

1~ = ~r* = "q ' P= — c~ = 1 " T ' ( 23 ) 

similar to those in (|14|) related to the semi-infinite slab. The Eqs. 123|) immediately 
follow from energy conservation considerations. 

A typical frequency dependence of finite photonic slab transmittance ()2rt[) is shown 
in Figs. H a, b, and c. For comparison, Fig. d shows the transmittance Ijl4|l of a 
semi-infinite photonic slab having the same periodic structure. The sharp peaks in 
transmittance in the vicinity of photonic band edge at larger N correspond to Fabry- 
Perot cavity resonances. At resonance, the electromagnetic field inside the slab is 
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D 



Figure 6. Light incident on a finite photonic slab of the thickness D. The 
subscripts /, -R, and P refer to the incident, reflected, and passed waves, 
respectively. The transmitted wave ^/t inside the slab may have Bloch componets 
propagating in either direction. 



close to a standing wave composed of one forward and one backward propagating Bloch 
eigenmodes with large and nearly equal amplitudes. The slab boundaries coincide with 
standing wave nodes, which determines the Bloch wavenumbers of the forward and 
backward components 

k s ^k g ±^s, 8 = 1,2,... (24) 

where k g corresponds to the photonic band edge. Usually, but not necessarily, k g 
equals or tt/L. The approximation (|24|) is valid if N 3> 1 and only applies to the 
resonances close enough to the photonic band edge. 

The dispersion function to (k) in the vicinity of a photonic band edge can be 
approximated as follows 



-jf (k-kg) 2 , where (25) 



2 



where 



The propagating mode group velocity u vanishes as lo — > uj g 

"^K^-^^k"^. (26) 

where ± corresponds to the forward and backward propagating waves, respectively. 
Inserting the resonance values (|24|l of the Bloch wavenumber into the dispersion 
relation Eq. (|25|l yields the resonance frequencies lu s as 

u^LOg-^-^-s) , .= 1,2,... (27) 
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Figure 7. Typical plots of transmittance versus frequency of lossless periodic 
stacks composed of different number N of unit cells. The frequency range shown 
includes a photonic band edge. The sharp transmission peaks near the band edge 
frequency are associated with Fabry-Perot cavity resonance. The case N = oo 
corresponds to a semi-infinite photonic slab and is similar to that shown in Fig. 
3. 



where uo g — uo (k g ) is the band edge. The dependence l|27|) is illustrated in Fig. [7]p. 

Let us focus on the Fabry-Perot cavity resonance closest to the photonic band 
edge. The respective frequency is 



^"*-f U • (28) 

At frequency u>i, the group velocities of the forward and backward propagating modes 
are 

7T ITUj" 

Ui«±u/'- = ± — 2-, (29) 

9 D LN y ' 

that is inversely proportional to the number N of the unit cells L in the slab. The 
resonance field amplitude inside the slab is proportional to the slab thickness 

* T (z) ps N^o sin (nzD- 1 ) + *i, < z < D, (30) 
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where an d ^1 are periodic functions of z comparable in magnitude with the incident 
wave ^ i. So, the maximum held amplitude is reached in the middle of the slab and 
is proportional to the slab thickness. The bandwidth Ai of the Fabry-Perot cavity 
resonance decreases sharply, as the number of unit cells increases 

^^t&^w (31) 

This is clearly seen in Fig. [7| 

At this point we would like to compare the frozen mode regime introduced in the 
previous subsection and the Fabry-Perot cavity resonance. Both effects result from 
coherent interference of light and can be thought of as photons trapped inside the 
periodic medium. Both effects are accompanied by a huge surge in electromagnetic 
held amplitude inside the photonic crystal. But that is where their similarity ends. 
Indeed, in the case of a Fabry-Perot cavity, the entire periodic stack (photonic slab) 
works as a resonator in which the trapped photons are spread all over the place. For 
this reason, all the major characteristics of Fabry-Perot cavity resonance are essen- 
tially dependent on the slab thickness. If the slab thickness D = NL is too large, then 
even small absorption or structural irregularity will completely smooth out the reso- 
nances. So, on the one hand, the slab should have enough layers to support distinct 
Fabry-Perot cavity resonances. But on the other hand, the number of layers should 
not be too large so that the losses and structural irregularities would not wipe out 
the effect. In addition, the number of layers essentially affects the resonance band- 
width. By contrast, the frozen mode regime is not a resonance in a usual sense of this 
word. Each trapped photon is now localized within certain small number of unit cells 
depending on the pulse bandwidth, while the slab size is not essential at all. Even if 
N — > oo, it does not affect any basic characteristics of the frozen mode regime, such 
as the bandwidth or the frozen mode amplitude. 



2.2.1. Photonic slab as a delay line In the case (|19[1 of a thin slab, the idea of a 
distinct pulse slowly propagating through the slab does not apply. On the other hand, 
one might be interested in the relation between the input and the output pulses, rather 
than in what is going on inside the photonic slab. This is the case, for example, if the 
photonic slab is used as a delay line. Let 'I'/ and fp be the input and output pulses, 
respectively, as illustrated in Fig. H3 The shape of the output pulse ^ p can be close 
to that of the input pulse ^fj regardless of whether or not the condition (|T$|l is met. 
If the shape of the pulse is indeed preserved, one can define the effective speed u of 
pulse propagation through the slab as 

fi = j, (32) 

where t is the transit time of the pulse passed through the slab. The quantity u is 
referred to as the group delay. The transit time determines the pulse delay due to the 
presence of the slab. Of course, in the case l(T§|) of a thick slab, the effective speed (l^2ll 
coincides with the pulse group velocity u. But now we consider the opposite situation 
H21fl . It turns out that under the resonance conditions, the transit time t of a thin 
photonic slab increases sharply, and the respective group delay u can be as low as 
10 -2 c, while the pulse passes through the slab with little reflection (see, for example, 
and references therein). In this sense, the pulse delay can be classified as a slow 



Slow light in photonic crystals 



16 



light effect, although the quantity u does not relate to the speed of any real pulse 
inside the photonic slab. In the rest of this section we briefly discuss this well known 
phenomenon. 

Let us estimate the group delay associated with Fabry-Perot cavity resonance. 
According to Eq. (|30|l . the electromagnetic energy Ti stored in the entire slab at the 
resonance is 

H oc |* T | 2 £> oc |*/| 2 LiV 3 . (33) 
This leads to the following rough estimate for the transit time t in (|32[1 

T/ 

ti oc — oc c^LN 3 . (34) 
bi 

The respective group delay (|3*2")l is 

Ul = j- oc — . (35) 

Note that if the number N is large, the value (|35|l of the group delay is much lower 
than the group velocity l|29[l of the propagating Bloch mode at the same frequency w\. 
The drawback, though, is that the bandwidth (|31|l of the Fabry-Perot cavity resonance 
shrinks even faster as the number N of unit cells increases. Eqs. (|35(l and (|31|l yield 
the following relation between the bandwidth Ai and the group delay u\ 

— oc — — . 36) 

u>i N c 

Comparison of the slow light bandwidth i|36|) with its ideal value (0} shows that the 
Fabry-Perot cavity resonance in a finite periodic photonic slab has a fundamental 
bandwidth disadvantage, if used as a delay line. 

Note that real optical delay lines are commonly based on periodic arrays of weakly 
coupled resonators, such as Fabry-Perot cavities, rather than on individual Fabry- 
Perot cavities (see, for example, |18l 1191 l2T)l 1211 l2~2l 1231 flQ and references therein). 

In conclusion, let us reiterate that in the cases other than (|19fl , there is no distinct 
pulse propagating inside the periodic medium and, therefore, the notion of slow light 
does not literally apply there. Further in this paper we assume that the condition (|19fl 
is satisfied, warranting the approximation of a semi-infinite photonic crystal. This 
allows us to investigate the slow light phenomenon in its pure form, when it is directly 
related to the speed of electromagnetic pulse propagation through the medium. In 
this case, the frozen mode regime associated with a stationary inflection point (1151) 
provides a unique possibility of converting a significant fraction of the incident light 
into a coherent mode with extremely low group velocity and drastically enhanced 
amplitude. 



3. Slow light in periodic layered media 

From now on we restrict ourselves to stratified media, which are periodic stacks of 
dielectric layers. Such systems are also referred to as photonic crystals with one- 
dimensional periodicity. A major reason for such a choice is that the electrodynamics 
of stratified media can be done within the framework of a rigorous analytical approach. 
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This is particularly important since the frozen mode regime involves a unique and 
spectacular behavior, so it would be desirable to be sure that such a behavior is not 
a numerical artifact. As soon as we assume that the semi- infinite photonic slab in 
Fig. |21 is a periodic array of plane-parallel uniform layers, we can give a much more 
detailed and meaningful description of the frozen mode regime. 

We start with some general remarks about electromagnetic eigenmodes in periodic 
layered media. Then we proceed to a semi-qualitative description of the situation 
taking place at different stationary points of the dispersion relations. A consistent 
and complete analysis based on the Maxwell equations will be presented in Sections 
5 through 12. 



3.1. Propagating and evanescent eigenmodes in periodic stacks of anisotropic layers 

Let \tj (z), (z) and ^>t {z) denote the incident, reflected and transmitted waves, 
respectively, as shown in Fig. [21 In the frequency domain, each of these waves can be 
explicitly represented by a column vector 



*(*) = 



E x {z) 

Ey{ Z ) 

H x {z) 

Hy(z) 



(37) 



where E x (z) , E y (z) , H x (z) , H y (z) are the transverse components of electromagnetic 
field. The exact definition of \& (z) is given in (|80|l and 1)81(1. The incident and reflected 
beams are plane monochromatic waves propagating in vacuum, while the transmitted 
electromagnetic field (z) inside the periodic layered medium is not a single Bloch 
eigenmode. At the slab boundary at z = 0, the three waves satisfy the standard 
boundary condition 

(0) + V R (0) = *t (0) , (38) 

implying continuity of the tangential field components 1(37(1 . Note that periodic stacks 
capable of supporting the frozen mode regime must include anisotropic layers with 
misaligned and/or oblique orientation of the principal axes. As a consequence, the 
reflected and transmitted waves in Fig. [3 will have an elliptic polarization even if the 
incident wave is linearly polarized. 

In the setting of Fig. [21 where the semi-infinite periodic layered array occupies 
the half-space z > 0, the transmitted wave (z) is a superposition of two Bloch 
components (Bloch eigenmodes) with different polarizations and different values of 
the Bloch wave number k. There are three possibilities. 

(i) Both Bloch components of the transmitted wave are propagating modes 

*r(«) = *pri(2)+^(2), z>0, (39) 

which means that the two respective values of k are real. For example, at 
uib < v < uj a in Fig. Q the transmitted wave is composed of two Bloch 
eigenmodes with two different real wave numbers k\ and Ai2 and two different 
group velocities u\ > and ui > 0. This constitutes the phenomenon of double 
refraction. 

(ii) Both Bloch components of are evanescent 

*T (z) = *e„l (z) + ^ev2 («) , Z > 0, (40) 



Slow light in photonic crystals 



18 



which implies that the two respective values of fc are complex with Imfc > 0. For 
example, this is the case when the frequency ui falls into the photonic band gap 
at u> > ujg in Fig. 2] The fact that Imfc > implies that the wave amplitude 
decays as the distance z from the semi-infinite slab surface increases. In the case 
(|40|l . the incident wave is totally reflected back to space by the semi-infinite slab, 
as seen in Fig. 03 

(iii) Of particular interest is the case where one of the Bloch components of the 
transmitted wave is a propagating mode with u > 0, while the other is an 
evanescent mode with Im fc > 

* T (z) = % r (z) + ^ev (Z) , Z>0. (41) 

For example, this is the case at the frequency range 

LU a < UJ < LOg (42) 

in Fig. ^ As the distance z from the slab/vacuum interface increases, the 
evanescent contribution ^ ev in (|41|l decays as exp (— zlmfc), and the resulting 
transmitted wave 'J't turns into a single propagating Bloch mode V? pr . 

Propagating modes with u < 0, as well as evanescent modes with Im fc < 0, never 
contribute to the transmitted wave inside the semi-infinite stack in Fig. [21 This 
fact is based on the following two assumptions: 

- The transmitted wave and the reflected wave are originated from the 
plane monochromatic wave incident on the semi-infinite photonic slab from 
the left, as shown in Fig. 

- The layered array in Fig. ^occupies the entire half-space and is perfectly periodic 
at z > 0. 

If either of the above conditions is violated, the field inside the periodic stack 
can be a superposition of four Bloch eigenmodes with either sign of the group velocity 
u of propagating contributions, or either sign of Imfc of evanescent contributions. For 
example, this would be the case if the periodic layered array in Fig. |2]had some kind 
of structural defects or a finite thickness like that presented in Fig. H3 

The propagating modes with u > and evanescent modes with Im fc > are 
referred to as the forward waves. Only forward modes contribute to the transmitted 
wave "J/t (z) in the case of a periodic semi-infinite stack. The propagating modes with 
u < and evanescent modes with Im fc < are referred to as backward waves. Since 
the backward Bloch waves are not excited in the setting of Fig. |3 they play no role 
in further consideration. 

In all three cases (|39l - l4T) | . the contribution of a particular Bloch eigenmode to 
the transmitted wave depends on the polarization "J/ of the incident wave. One 
can always choose the incident wave polarization so that only one Bloch component 
is excited. In such a case, is a single Bloch eigenmode. 

Only propagating modes contribute to the normal component St of the energy 
flux inside a periodic semi-infinite slab. Evanescent modes do not participate in energy 
transfer in such a case. In the important particular case of a single propagating mode 
(* T = * pr ), we have from (JT3J) and ((HI) 



St = Wu = tSi, 



(43) 
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where W ~ i s the energy density associated with the transmitted propagating 

mode. 

The assumption that the transmitted wave (z) is a superposition of 
propagating and/or evanescent Bloch eigenmodes may not be valid at stationary points 
11211 of electromagnetic dispersion function u(k), because each stationary point is a 
degeneracy point of the frequency spectrum. For example, if the frequency u> exactly 
coincides with stationary inflection point defined by (|15|) . the transmitted wave H>t (z) 
is dominated by a (non-Bloch) Floquet eigenmode linearly growing with z, which 
constitutes the frozen mode regime j^HEOl- At all other frequencies, the transmitted 
wave is a superposition of two forward Bloch modes, each of which can be either 
propagating or evanescent. A detailed analysis of this and related phenomena is 
presented further in this paper. 

Knowing the eigenmode composition of the transmitted wave (z) we can give 
a semi-qualitative description of what happens when the frequency u> of the incident 
wave approaches one of the stationary points (|12fl in Fig. ^ A complete analysis 
based on the Maxwell equations will be presented later in the paper. 



3.2. Photonic band edge 

We start with the simplest case of a photonic band edge. Just below the band edge 
frequency uj g in Fig. ^ the transmitted field (z) is a superposition (|41|) of one 
propagating and one evanescent Bloch component. Due to the boundary condition 
H38fl at the slab/vacuum interface, the amplitude of the transmitted wave at z = is 
comparable with that of the incident wave. In the case of a generic polarization of the 
incident light, the amplitudes of the propagating and evanescent Bloch components 
at z — are also comparable to each other and to the amplitude of the incident light 

IV (0)| ~ \Vev (0)| ~ |*/| , at u < u g . (44) 

As the distance z from the slab surface increases, the evanescent component ^ ev (z) 
decays rapidly, while the amplitude of the propagating component remains constant. 
Eventually, at a certain distance from the slab surface, the transmitted wave (z) 
becomes very close to its propagating Bloch component 

f T (z) S3 V pr (z), at z > L, oj < uj g . (45) 

The evanescent component $ ev of the transmitted wave does not display any 
singularity in the vicinity of u> g . By contrast, the propagating mode ^ pr develops a 
singularity as u> — * u> g — 0, which is associated with vanishing group velocity as in Ijl2|l . 
At ui > u g , the propagating mode turns into another evanescent mode in 1)40(1 . 

The dispersion relation in the vicinity of the band edge g in Fig. ^ can be 
approximated as 

u>" 

Ug ~ CJ S3 -± (k - kg) , U^Ulg. 

This yields the following frequency dependence of the propagating mode group velocity 
u below the photonic band edge 

U=^^^{kg-k)^^ j {u J g-Uj) l/ \ U<UJg. (46) 
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The energy flux l|43|> associated with the slow propagating mode is 

S T & ( W V^( UJ 9^ UJ ) 1/2 ^ at _ ( 47 ) 

\ 0, at lo > uj g 

where 

The latter estimation follows from (|44|l and applies to the case of a generic polarization 
of the incident wave. The semi-infinite slab transmittance (|f4|) in the vicinity of ui g is 



S/ I 0, at w > w 9 

where according to Ij48(l 

W |*/| 2 



(49) 



Si Si 

The relation l|49|) is illustrated by the numerical example in Fig. [3J 

Equation (|49|l expresses the well-known fact that in the vicinity of an 
electromagnetic band edge, the semi-infinite photonic crystal becomes totally 
reflective, as illustrated in Fig. [3] This implies that as lo — > Lo g , only an infinitesimal 
fraction of the incident light energy is converted into the slow mode. 



3. 3. Other extreme points of spectral branches 

For specificity, let us consider the stationary point a of the dispersion relation in Fig. 
n which qualitatively is not different from the point b. At frequencies right below 
ui a , the transmitted wave is a superposition l|39|) of two propagating eigenmodes, 
one of which is the slow mode and the other is a regular forward propagating mode. 
The slow mode develops a singularity at lo — to a similar to that of the respective 
slow mode in the vicinity of the band edge frequency u) g , while the other propagating 
mode (the fast mode) remains regular in the vicinity of Lo a and does not produce 
any anomaly. The two forward modes contribute additively to the energy flux St, 
but the contribution of the fast mode remains regular in the vicinity of u) a , while the 
contribution of the slow mode shows the same singular behavior as that described by 
Eqs. (|4"7| and jl5jl. Fig. Olprovides a graphic illustration of such a behavior. 

The important point is that similar to the situation in the vicinity of a photonic 
band edge, at to = ui a and lo = u>b the contribution of the respective slow mode to the 
transmitted wave vanishes. In other words, in terms of slow mode excitation, the 
stationary points a and b in Fig. ^are no different from the band edge g. 



3.4- Stationary inflection point: the frozen mode regime 

A sharply different situation develops in the vicinity of a stationary inflection point 
11511 of the dispersion relation (point in Fig. 2J. According to l|15|) . the dispersion 
relation in the vicinity of loq can be approximated as follows 

lo'" 

co-io^-^ik-kof, (50) 
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where 

j„ = ( d 3 iv 

, <9fc 3 / k=k 

The propagating mode group velocity u vanishes as lu approaches too 

® U ~ 1 W f\ i ^2 _ 6 2/3 , W xl/3 / n2/3 

Sfc 2 ~ — (w-^o) ■ (51) 

But remarkably, the electromagnetic energy density W associated with the transmitted 
frozen mode diverges as ui — > ujq 

tj/ 2r5/ . ,,,.-1/3 , \-2/3 /r-r,\ 

W ~ "^273" ( w o ) ' - cj ) , (52) 

where Sj is the fixed energy flux of the incident wave. The slab transmittance r 
remains finite even at w = wo, as illustrated in Fig. |3] As a result, the energy flux 
l)43|l associated with the transmitted frozen mode also remains finite and can even 
be close to unity in the vicinity of luq. The latter implies that the incident light is 
completely converted to the frozen mode with infinitesimal group velocity Q51JI and 
diverging energy density ifS^l) . 

Let us consider the structure of the frozen mode. At w k wo, the transmitted 
wave is a superposition i|41|) of one propagating and one evanescent Bloch 
component. In contrast to the case of a photonic band edge, in the vicinity of ujq both 
Bloch components of *S?t develop strong singularity. Specifically, as the frequency ui 
approaches wo, both contributions grow sharply, while remaining nearly equal and 
opposite in sign at the slab boundary 

^ pr (0) w -^ ev (0) oc \u> - cj r 1/3 , as w w . (53) 

Due to the destructive interference 1)53(1 . the resulting field 

*T (0) = «V (0) + *«, (0) 

at the slab boundary is small enough to satisfy the boundary condition (|38JI . as 
illustrated in Fig. 03 As the distance z from the slab boundary increases, the 
evanescent component ^ ev decays exponentially 

* el) (z) « y ev (0) exp (-z Im k) 

while the amplitude of the propagating component ^ pr remains constant and very 
large. As a consequence, the amplitude of the resulting transmitted wave \1/t (z) 
sharply increases with the distance z from the slab boundary and, eventually, reaches 
its large saturation value corresponding to the propagating component \t pr , as 
illustrated in Fig. El 

If the frequency to of incident light is exactly equal to the frozen mode frequency 
w , the transmitted wave *S?t (z) does not reduce to the sum (|4T|l of propagating 
and evanescent contributions 1211130]. Instead, it corresponds to a non-Bloch Floquet 
eigenmode diverging linearly with z 



tSt 

* T (z) - * T (0) cx z J — * , atw = w . 
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Figure 8. Destructive interference of the propagating and evanescent 
contributions to the resulting field at the slab/ vacuum interface under 

the frozen mode regime: a) resulting field amplitude Iv&y (0)| 2 , b) amplitude 
I *&pr (0) | 2 of the propagating component, c) amplitude \9 e v (0) | of the evanescent 
component. 



a) Resulting field b) Propagating component c) Evanescent component 
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Figure 9. The transmitted electromagnetic field 1411 and its propagating and 
evanescent components inside semi-infinite slab in close proximity of the frozen 
mode regime: (a) the amplitude \^>t (z)\ 2 of the resulting field, (b) the amplitude 
|\E , pr (2)| 2 of the propagating contribution, (c) the amplitude \9 ev (z)\ 2 of the 
evanescent contribution. Due to destructive interference of the propagating and 
evanescent components, the resulting field amplitude at z = is small enough to 
satisfy the boundary conditions (38). The amplitude |\Pj| 2 of the incident wave 
is unity. The distance z from the slab boundary is expressed in units of L. 
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Evidently, the frozen mode regime associated with stationary inflection point (|15fl 
provides an ideal and unique situation in terms of slow mode excitation. Indeed, in 
this case virtually all the incident light energy can be converted into the slow mode 
with greatly enhanced amplitude. 

A consistent mathematical analysis of the frozen mode regime is rather 
sophisticated and will take a great deal of our attention further in this paper. 
Specifically, the fundamental fact that at ui — uj$, the energy flux of the frozen mode 
remains finite in spite of vanishing group velocity, is rigorously proven in Section 9 
(see Eqs. <|421[) through l|449|) and related explanations). 

3.5. Degenerate band edge 

While the situation with the regular band edge appears quite obvious and not 
particularly interesting from the perspective of slow light, the so-called degenerate 
band edge proves to be quite different |34| . An example of an electromagnetic 
dispersion relation with degenerate band edge is shown in Fig. 



a) Regular band edge b) Degenerate band edge c) Double band edge 




123456 123456 123456 

Wave number k Wave number k Wave number k 



Figure 10. Dispersion relations lu (k) of three periodic stacks with slightly 
different layer thicknesses. The plot in the middle displays a degenerate band 
edge (54). 

In the vicinity of the degenerate band edge d, the dispersion relation can be 
approximated as 

u d - u « ^ (k - k d f , where W »" = ^ ^ . (54) 

Similarly to the regular band edge, below the degenerate band edge frequency u>d, the 
transmitted field \I/t is a superposition (|41|) of one propagating and one evanescent 
components, while above ujd, the transmitted wave is a combination (|40|) of two 
evanescent components. The critical difference though is that now both the Bloch 
components display strong singularity. Specifically, as the frequency ui approaches w^, 
both Bloch contributions grow sharply, while remaining nearly equal and opposite in 
sign at the slab boundary 

*pr (0) « -V e v (0) oc \uj d ~ u\~ 1/4 , as w ^ w,, - 0. (55) 
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The destructive interference l|55|l ensures that the boundary condition i|38|) is satisfied, 
while both Bloch contributions to "J^ (z) diverge. As the distance z from the slab 
boundary increases, the evanescent component ty ev (z) dies out, while the propagating 
component ty pr (z) remains huge. Eventually, at z 3> L, the resulting transmitted 
wave *5>t (z) coincides with the propagating Bloch eigenmode ^ pr (z). If the frequency 
u of incident light is exactly equal to u>d, the transmitted wave (z) does not reduce 
to the sum (|41|) of propagating and evanescent contributions. Instead, it corresponds 
to a non-Bloch Floquet eigenmode diverging linearly with z 

<f T (z) - * T (0) oc at u = uj d . 

The above behavior appears to be very similar to that of the frozen mode regime 
described in the previous subsection. In both cases, Figs. 03 and provide a good 
graphical illustration of the electromagnetic field distribution inside the slab in the 
vicinity of the relevant stationary point. Yet, there is a crucial difference between the 
two cases. In spite of its huge diverging amplitude H55[). the transmitted wave ^>t 
does not provide any energy flux in the immediate proximity of a degenerate band 
edge. Indeed, according to l|55ll . as u> approaches ujd, the energy density W of the 
transmitted wave diverges as 

W OC l^prl 2 \u>d — U>\ , at U) LL>d- 

But, from Eq. I|54(l one can derive that the respective slow mode group velocity 
vanishes even faster 

w 'd" n i \3 24 3 / 4 1/4 3 / 4 
u ~ -f- \ k - k d) ~ — — ) - w ) at w $ w,j 

6 6 

As the result, the energy flux of the transmitted wave vanishes, as one approaches the 
degenerate band edge 

a iir f (^d — ^) 1 ^ 4 , at U> ^ bJd 

^ U, at u> > ujd 

and so does the slab transmittance r 

t ^ f (u d - uj) 1/4 , at w £ uj d 
\ 0, at w > uod 

By contrast, in the case of the frozen mode regime the slab transmittance remains 
finite and a significant fraction of the incident light energy goes to the slow mode. 

The situation at a degenerate band edge can be viewed as intermediate between 
the frozen mode regime and the vicinity of a regular band edge|34j. Indeed, on the one 
hand, the incident wave at to = ujd is totally reflected back to space, as would be the 
case at a regular band edge. On the other hand, the transmitted field amplitude inside 
the slab becomes huge as u — > uj d , which is similar to what occurs in the frozen mode 
regime. The large amplitude of the transmitted wave at u) ~ u>d can be very attractive 
for a variety of practical applications, although such a behavior cannot be qualified 
as a slow light case. Detailed analysis of some peculiar electromagnetic properties 
associated with degenerate frequency band edge (|54|l can be found in |331 13~4*| . 
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4. Physical conditions for the frozen mode regime in layered media 

The frozen mode regime is associated with a stationary inflection point ljl5|l of 
the electromagnetic dispersion relation. Leaving the proof of this statement to the 
following sections, here we establish the conditions under which the dispersion relation 
of a periodic layered array can develop the singularity (JTSJ. We will see that only 
special layered structures incorporating anisotropic layers can display this property. 
In the following sections, based on the Maxwell equations, we will show that indeed 
the stationary inflection point (|15|l is uniquely associated with the frozen mode regime. 

4-1- Axial dispersion relation: basic definitions 

We start with the generalization of the frozen mode concept to the case of oblique 
light incidence. 

Consider a monochromatic plane wave obliquely incident on a periodic semi- 
infinite stack, as shown in Fig. [3 Let 'J'/, and *S?t denote the incident, reflected 
and transmitted waves, respectively. Due to the boundary conditions l|38|) . all three 
waves Vf/, and 'J't must be assigned the same pair of tangential components 
k x , k y of the respective wave vector [2] 




while their axial (normal) components k z can be different. Hereinafter, the normal 
component of the transmitted Bloch waves propagating inside the periodic layered 
medium will be referred to as the wave number and denoted by the symbol fc, rather 
than k z , so that inside the periodic stack (at z > 0) 

k = {k Xl ky,k) . (57) 

Unlike k x and k y , the z component k of the Bloch wave vector (|57|) is defined up to a 
multiple of 2tt/L 

k = k + 2wN/L, (58) 

where L is the period of the layered structure and N is an integer. For given k x , k v 
and to, the value k is found by solving the time-harmonic Maxwell equations l|79|) in 
the periodic medium, as will be done in the following sections. The result can be 
represented as the axial dispersion relation, which gives the relation between lu and k 
at fixed k x , k y 

lu = lu (k) , at fixed k x , k y . (59) 

It can be more convenient to define the axial dispersion relation as the relation between 
uj and k at fixed direction n of incident light propagation 

lu = lu (k) , at fixed n x , n y , (60) 

where the unit vector n can be expressed in terms of the tangential components (|56|l 
of the wave vector 



n x = k x c/u>, n y = kyc/uj, n z = J 1 - (n 2 x + n y ). (61) 
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4-1.1. Axial stationary inflection point and the frozen mode regime Suppose that 
at k — ko and u> = ojq = lo , one of the axial spectral branches l|59|l develops a 
stationary inflection point for given (k x , k y ), namely 



The value 



(63) 



in Eq. 1(62(1 is the axial component of the group velocity, which vanishes at k = kg. 
Observe that 

fr) ™ d "» = (^) • ,64 » 

^ UK x/k,k y \ UK V/k,k x 

representing the tangential components of the group velocity, may not be zero at 
k — fco. The spectral singularity (|62|l is called the axial stationary inflection point. 

One can also use another definition of axial stationary inflection point (|62|l . which 
is based on the axial dispersion relation (|60|) rather than l|59|) . namely 

at« = «o and n = no: ^'(^J = °> ( J * °' 

(65) 

The partial derivatives in l|65(l are taken at constant (n x ,n y ), rather than at constant 
(k x ,k y ). Both definitions 1(62(1 and 1(65(1 are equivalent to one other. In the particular 
case of normal incidence in which n || k \\ z, the axial stationary inflection point (|62|l 
or, equivalently, l|65|l turns into a regular stationary inflection point l|15(l . 

The axial frozen mode regime associated with the singularity 1)62(1 is very similar 
to its particular case, the regular frozen mode regime, related to the regular stationary 
inflection point 115(1 . Specifically, in the axial frozen mode regime, obliquely incident 
light can enter the semi-infinite photonic crystal with little reflection, where it is 
completely converted into a coherent mode with infinitesimal normal component ((63(1 
of the group velocity and drastically enhanced amplitude. The energy density of the 
axial frozen mode displays the same resonance-like behavior 1(52(1 . The only difference 
between the axial and the regular frozen mode regime is that in the former case, 
the tangential component ((64|) of the group velocity remains finite at u> — uiq. The 
specificity of the axial frozen mode regime as compared to the regular one is discussed 
in |29j . Further in this paper we will focus exclusively on the common features of 
these two cases. Either of them will be referred to simply as the frozen mode regime. 



4-2. Spectral asymmetry in periodic stacks 

The (axial) stationary inflection point is indeed associated with the frozen mode 
regime. But not every periodic layered media can display such a spectral singularity. 
It turns out that a necessary condition for the existence of an axial stationary inflection 
point and, therefore, a necessary condition for the frozen mode regime is the following 
property of the electromagnetic dispersion relation of the periodic stack 

lo (k x ,k y ,k) ^ lo (k x ,k y ,-k) or, equivalently, uo (n x ,n y ,k) ^ u (n x ,n y ,-k) . (66) 
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The property l|66|) is referred to as axial spectral asymmetry. Further in this 
paper, we will use the simplified notation l|59|) for the axial dispersion relation. In this 
notation, the requirement 1|66[1 of axial spectral asymmetry takes the following form 

lj (k) ? u (-k) , (67) 

where k is the z component (|5H|) of the Bloch wave vector k. A robust frozen mode 
regime only occurs if the degree of spectral asymmetry (|67|) is significant. For brevity, 
hereinafter, the quantity k will be referred to as the Bloch wave number, although in 
the case of oblique propagation, k is just the normal component of the Bloch wave 
vector k. 

In the particular case of normal wave propagation, in which ft \\ k \\ z, the 
requirement Ij66|l of axial spectral asymmetry reduces to 

w (jcj ^ oj (-fc) ,k\\ z. (68) 

This kind of asymmetric dispersion relation can occur only in periodic structures with 
some of the constitutive components being magnetic and displaying nonreciprocal 
Faraday rotation Significant spectral asymmetry requires strong Faraday 

rotation. The simplest periodic array supporting the spectral asymmetry ()68|l is shown 
in Fig. ITT1 



1 2 F 




L 



Figure 11. The simplest periodic stack supporting asymmetric dispersion 
relation (68). A unit cell L of this stack comprises three layers: two anisotropic 
layers 1 and 2 with misaligned in-plane anisotropy (the A - layers), and one 
magnetic layer F with magnetization shown by the arrows. 



At microwave frequencies, there exist a number of magnetic materials displaying 
low losses and strong Faraday rotation. But, at infrared and optical frequencies, the 
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magnetic materials with sufficiently strong Faraday rotation are usually too lossy for 
our purposes. Therefore, if we are interested in optical frequencies, we have to rely on 
non-magnetic stacks, in which the regular spectral asymmetry (|68|l is impossible. By 
contrast, the axial spectral asymmetry (|66|l or (|67|) does not require magnetic layers 
and can occur even in perfectly reciprocal non-magnetic stacks. 

The physical conditions under which the electromagnetic dispersion relation of a 
non-magnetic layered structure can develop the (axial) spectral asymmetry (|67[) and, 
thereby, support the (axial) frozen mode regime can be grouped in two categories. The 
first one comprises several symmetry restrictions. The second category includes some 
basic qualitative recommendations which would ensure the robustness of the frozen 
mode regime, provided that the symmetry conditions for the regime are met. In what 
follows we briefly describe those conditions and then show how they apply to periodic 
stacks incorporating some real dielectric materials. 

There are two fundamental necessary conditions for the frozen mode regime in 
a non-magnetic (reciprocal) periodic layered structure. The first one is that the 

Bloch dispersion function lu (jtj in the periodic layered medium must display the 

axial spectral asymmetry (|67|) . This condition is necessary for the existence of the 
axial stationary inflection point (|62fl in the electromagnetic dispersion relation of 
an arbitrary periodic layered medium. The second necessary condition is that for 
the given direction k of wave propagation, the Bloch eigenmodes with different 
polarizations must have the same symmetry. In the case of oblique propagation in 
periodic layered media, the latter condition implies that for the given k, the Bloch 
eigenmodes are neither TE nor TM: 

is neither TE nor TM. (69) 

The condition l|t)7[) imposes certain restrictions on (i) the point symmetry group 
G of the periodic layered array and (ii) on the direction k of the transmitted wave 
propagation inside the layered medium, while the condition (|69|) may impose an 
additional restriction on the direction of k. 

The restriction on the symmetry of the periodic stack stemming from the 
requirement l|67|) of the axial spectral asymmetry is 

m z £G and 2 Z (£ G. (70) 

where m z is the mirror plane parallel to the layers, 2 Z is the 2-fold rotation about 
the z axis. An immediate consequence of the criterion l|7U|) is that at least one of the 
alternating layers of the periodic stack must be an anisotropic dielectric with 

e xz ^ and/or e yz ^ (71) 

where the z direction is normal to the layers. Otherwise, the operation 2 Z will be 
present in the symmetry group G of the periodic stack. 

The simplest and the most practical example of a non-magnetic periodic stack 
satisfying the criterion H7(J[) is shown in Fig. It is made up of anisotropic A layers 
alternating with isotropic B layers. The respective dielectric permittivity tensors are 



Eyy 

&xz &ZZ 



£b = 



£ B 
e B 
e B 



(72) 
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Figure 12. Periodic layered structure with two layers A and B in a unit cell 
L. The A layers (hatched) are anisotropic with one of the principle axes of the 
dielectric permittivity tensor e making an oblique angle with the normal z to the 
layers (e xz ^ 0). The B layers can be isotropic. The x — z plane coincides with 
the mirror plane m y of the stack. 

For simplicity, we assume 

P-A = As = I- ( 73 ) 
The stack in Fig. 1121 has the monoclinic symmetry 

2 y /m y (74) 

with the mirror plane m y normal to the y - axis. Such a symmetry is compatible with 
the necessary condition (|70|l for the axial spectral asymmetry (|67|) . Therefore, the 
periodic array in Fig. 1121 can support the frozen mode regime, provided that incident 
beam lies neither in the x — z, nor in y — x plane |29j 

n x 7^ and n y =/= 0. 

If all the above necessary conditions are met, then the (axial) frozen mode regime 
is, at least, not forbidden by symmetry. More details on the symmetry aspects of the 
frozen mode regime can be found in jSU], Section II and Sections I and II. 

In practice, as soon as the symmetry conditions are met, one can almost certainly 
achieve the (axial) frozen mode regime at any desirable frequency lo within a certain 
frequency range. The frequency range is determined by the layer thicknesses and the 
dielectric materials used, while a specific value of lo within the range can be selected 
by the direction ft of the light incidence. The problem is that unless the physical 
parameters of the stack layers lie within a certain range, the effects associated with 
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the frozen mode regime can be insignificant or even practically undetectable. The 
basic guiding principle in choosing appropriate layer materials are discussed in Ref. 

H3EI] 

The biggest challenge at optical frequencies lies in the fact that most of the 
commercially available optical anisotropic crystals have weak anisotropy. According 
to PHI, this would push the axial stationary inflection point (|62|) very close to the 
photonic band edge and make the photonic crystal almost 100% reflective. This indeed 
would be the case if we tried to realize the frozen mode regime at the lowest frequency 
band. But, in Ref. |3T] it was shown that the above problem can be successfully 
solved by moving to a higher frequency band. So, a robust axially frozen mode 
regime with almost complete conversion of the incident light into the frozen mode 
can be achieved with the commercially available anisotropic dielectric materials such 
as NbLi0 3 , YV0 4 , etc. 

5. Electrodynamics of lossless stratified media 

This section starts with a description of some basic electrodynamic properties of 
stratified media composed of lossless anisotropic layers. Then we turn to the important 
particular case of unbounded periodic layered arrays (periodic stacks), where the 
electromagnetic eigenmodes are Bloch waves. Then we consider the problem of 
electromagnetic energy flux in lossless stratified media. Finally, we outline the 
electromagnetic scattering problem for a semi-infinite photonic slab. The material 
presented in this section is sufficient for a numerical analysis of slow light phenomena 
in periodic layered media. This approach was used in |291 I3UI 131) to analyze 
the frozen mode regime in magnetic and non-magnetic periodic stacks. Yet, to 
develop a consistent analytical picture of the frozen mode regime we shall need a 
more sophisticated mathematical framework based on a perturbation theory for non- 
diagonalizable degenerate matrices. This problem will be addressed in the following 
sections. 

5.1. Reduced time-harmonic Maxwell equations 

Our treatment is based on the time-harmonic Maxwell equations in heterogeneous 
nonconducting media 

V x E(r) =i-B(f), V x H(r) = -i-TD (r) . (75) 
c c 

Electric and magnetic fields and inductions in Eq. (|75|l are related through the linear 
constitutive equations 

D (r) = g (r) E (r) , B (r) = £ (r) H (r) . (76) 

All variables in Eqs. Q75f) and l|76ll are frequency dependent. In lossless media, the 
material tensors e (f) and fi (r) are Hermitian 

In lossless media: e (r) = (r) , fi(r) = (J (r) , (77) 

where the dagger f denotes the Hermitian conjugate. In lossless non-magnetic media, 
both tensors e and jx are also real and symmetric 

In lossless non-magnetic media: e (r) = e* (r) = e T (f) , fx(r) = f 1 * (r) = (i T (r) , 

(78) 
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where the asterisk denotes the complex conjugate and the superscript T denotes matrix 
transposition. In magnetically polarized lossless media, the Hermitian material tensors 
<|77[1 may have a skew-symmetric imaginary part which is responsible for the non- 
reciprocal effect of Faraday rotation [2] . 

In a stratified medium, the second rank tensors e (r) and (i (f) depend on a single 
Cartesian coordinate z, and the Maxwell equations (|75|l can be recast as 



V x E (r) = i-fi (z) H (f) , V x H (r) = -i-e (z) E (f) , 
c c 

Solutions for Eq. I|79ll can be sought in the following form 

E (r) = e *(k*x+k y y)E ( z ) , H (r) = e i( - k " x+k ^H (z) , 



(79) 



(80) 



which can be interpreted as the "tangential" Bloch representation. The substitution 
H80JI allows separation of the tangential components of the fields into a closed system 
of four linear ordinary differential equations 



d z ^! (z) = i-M (z) * (z) , ^ (z) 



E x {z) 

Ey(z) 

H x {z) 

Hy(z) 



(81) 



where the 4x4 matrix M (z) is referred to as the (reduced) Maxwell operator. The 
normal held components E z and H z do not enter the reduced Maxwell equations (|81ll 
and can be expressed in terms of the tangential field components from Eq. 18111 as 



E z ( n x Hy -\- HyH x £ XZ E X Sy Z Ey^j & zz , 
H z = (n x E y - n y E x - fi* xz H x - ^* yz H y ) [i zz , 

where n x — ck x /uj, n y — ck y /uj. 

The explicit expression for the Maxwell operator M (z) in Eq. H81(l is 



(82) 



M(z) 



M 2 i M 22 



(83) 



where 
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In the important particular case of k x 
simpler form 

M{z) 



k y = 0, the Maxwell operator l|83|l has a 





M 21 







(84) 



Only in this case of k || z, the fields E (z) and H (z) coincide with the actual electric 
and magnetic fields E (r) and H (r) , as follows from the relation H80|) . 
The (reduced) Maxwell operator M (z) is a function of: 

- the local values of material tensors e (z) and fi (z), 

- the space coordinate z, via the material tensors e (z) and fj, (z), 

- the tangential components k x , k y of the wave vector, 

- the frequency lu. 

Different versions of the reduced Maxwell equation l|81|l can be found in the 
extensive literature on electrodynamics of stratified media (see, for example, and 
references therein). 

The 4x4 matrix M (z) in Eq. (|83|l has the fundamental property of J-Hermitivity 
defined as 

M f = JMJ, 



where 



J = J f = J~ 






-1 

1 





-1 





(85) 



(86) 



The property 1)85(1 relates to conservation of electromagnetic energy in lossless media, 
it plays an important role in further consideration. 



5. 2. The transfer matrix formalism 
The Cauchy problem 

d z V (z) = i-M (z) * (z) , * (z ) = * (87) 
c 

for the reduced Maxwell equation l|81|l has a unique solution 

9(z) =T(z,z )*(z ) (88) 

where the 4x4 matrix T (z, zq) is referred to as the transfer matrix. From the definition 
(HBJ it follows that 

T (z, z ) = T (z, z') T (z', z ) , T (z, z ) = T' 1 (z , z) , T (z, z) = /, (89) 

where / is the identity matrix. The transfer matrix T (z, zq) allows determination of 
the time-harmonic electromagnetic field ^ (z) at an arbitrary point z of the stratified 
medium once its value ^ (zq) at any particular point zq. 

The matrix T (z, zo) itself is uniquely defined by the following Cauchy problem 

d z T (z, z ) = i-M (z) T (z, z ) , T (z, z) = I. (90) 
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The equation l|9l)|l . together with J-Hermitivity (|85|) of the Maxwell operator M (z), 
implies that the transfer matrix T (z, z ) is J -unitary 

T t = j T -i J; ( 91 ) 

as shown in Ref. [3D] . The J-unitarity (|91|l of the transfer matrix T = T(z,zo) 
imposes the following constraint on its set of four eigenvalues Ci, i = 1, 2, 3, 4 

id , e 2 > cs , en = {cr 1 , G" 1 > C3" 1 , a 1 } , (92) 

which also implies that 

|detT| = l. (93) 

5.2.1. The transfer matrix of a stack of uniform layers The greatest advantage of 
the transfer matrix approach stems from the fact that the transfer matrix T$ of an 
arbitrary stack of layers is a sequential product of the transfer matrices T m of the 
constituent layers 

T s = l[T m . (94) 

m 

According to Eq. I|9U|I . if each individual layer m is homogeneous, the corresponding 
single-layer transfer matrices T m can be explicitly expressed in terms of the respective 
Maxwell operators M m 



exp (l^-ZmM^ , (95) 



T,,. 

c 



where z m is the thickness of the m-th layer. The explicit expression for the Maxwell 
operator M m of an arbitrary uniform layer of anisotropic dielectric material is given 
by Eq. |(H3J|. Thus, Eq. together with (|93|) and (|55)l . give an explicit expression 

for the transfer matrix T$ of an arbitrary stack of anisotropic dielectric layers. The 
elements of the J-unitary matrix T$ are functions of: 

- the material tensors e and fx in each layer of the stack, 

- the layer thicknesses d m , 

- the frequency U), 

- the tangential components k x = —n x , k y = ^n y of the wave vector. 
In the case of k x — k y — 0, we have instead of Eq. 193(1 

if k x = k y = 0, det T s = 1, (96) 

that can be derived directly from Eqs. I|95() and l|83|l . 

5.3. Periodic layered arrays. Bloch eigenmodes. 

In a periodic layered medium, all material tensors are periodic functions of z 

e(z) = e (z + L) , ft, (z) = p, (z + L) , 
and so is the Maxwell operator M(z) in Eq. I|81|) . 

M{z + L)= M [z) , (97) 
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where L is the length of a unit cell of the periodic stack. A Bloch solution (z) of 
the reduced Maxwell equation (|81|l with the periodic operator M(z) should satisfy the 
following relation 

* fe (z + L) = e lkL V k (z) , (98) 
where k is the normal component of the Bloch wave vector 

k = k z . (99) 

Unlike k x and k y , the z component l(99() of the Bloch wave vector is defined up to a 
multiple of 2n/L 

k = k + 2ttN/L, 

where N is an integer. Hereinafter, the normal component k z of the Bloch wave vector 
k will be referred to simply as the wave number and denoted with symbol k, rather 
than k z . 

The definition of the T-matrix together with Eq. give 

* fc (z + L) = T (z + L, z) * fe (z) = e lkL V k (z) . (100) 
Introducing the transfer matrix of a primitive cell 

T L =T(L,0) (101) 

we have from Eq. Ij 100(1 

T L $ fe = e 2fei $ fc , where $ fc = (0) . (102) 

Thus, the eigenvectors of the transfer matrix of the unit cell are uniquely related 
to the eigenmodes {z) of the reduced Maxwell equation (|81|l through the relations 

$ i = *i(0), *= 1,2,3,4. (103) 

The respective four eigenvalues 

C^e^, i = 1,2,3,4 (104) 
of Tl are the roots of the characteristic equation 

F (C) = 0, where F (Q = det (t l - (i) . (105) 



For ciny given to Rnd {kxi 

), the characteristic equation l|105|) defines a set of four 
values C 2 i C3 7 d}, or equivalently, {fci, k%, k%, ^4}. Real k (or, equivalently, |£| = 1) 
correspond to propagating Bloch waves (propagating modes), while complex k (or, 
equivalently, |£| 7^ 1) correspond to evanescent modes. Evanescent modes are relevant 
near photonic crystal boundaries and other structural irregularities. 

The J-unitarity (|91|l of Tl imposes the restriction ((92JI on the eigenvalues ((104(1 . 
which can be recast as 

{h ,k 2 ,k 3 ,h} = {k{ , k* ,k* 3 ,kt}, (106) 

for any given w and k x ,k y . In view of the relations l(106|) or l(92|) . one can distinguish 
the following three different situations. 
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(i) All four wave numbers are real 

ki = k\, k 2 = k2, k 3 = fcg, fot = k\. (107) 

or, equivalently, 

lCil = IC 2 | = IC 3 | = IC4l = i- (108) 

The respective four Bloch eigenmodes are propagating. 

(ii) Two real and to complex wave numbers 

ki — k^, k2 — k* 2 i ki = fcg, where k% ^ k^, k^^k\. (109) 
or, equivalently, 

Kil = |Cal = 1, U = l/a where |C 3 | , IC 4 I * 1- (HO) 

Two of the four Bloch eigenmodes are propagating and the remaining two are 
evanescent with complex conjugate wave numbers. 

(iii) All four wave numbers are complex 

^2 — k^ = fcg, where k\ ^ k\, k 2 ^ k 2 , k 3 ^ fcg, k^^k\. (HI) 
or, equivalently, 

C 2 = 1/Ci, C 4 = 1/CI, where |Cxl . Kal , Ksl . IC 4 I * !• (112) 

This situation relates to a frequency gap, where for given u> and k x ,k y , all four 
Bloch eigenmodes are evanescent. 

5.4- Symmetry of the dispersion relation 

Below we will see that the dispersion relation is symmetric under k — > — k if and only 
if the transfer matrix Tl is similar to its inverse. Indeed, assume that 

T L = V- x Tl x V (113) 

where U is a nonsingular 4x4 matrix. This assumption together with the property 
H91fl of J-unitarity, imply also the similarity of Tj, and Tj 

T L = V^TlV, (114) 

where V = JU. Either of the above two relations imposes the following additional 
restriction on the eigenvalues (|l(J4|t of Tl for given uj and k x , k y 

{Ci.Ca.Cs.CJ = {Cl\G\G\Q 1 h (H5) 

or, equivalently, 

{^1,^2,^3,^4} = {-ki,-k 2 ,-k 3 ,-k4}- (116) 

The relation (|1 lfi|> is referred to as axial spectral symmetry. It applies both to 
propagating and evanescent solutions. 

Let us consider the symmetry relation l|lltj|) in more detail. Assume that k\ = fc* 
is a real wave number corresponding to a propagating eigenmode. The relation l|116fl 



Slow light in photonic crystals 



36 



implies that for given oj and k x ,k y , there is another real wave number fc 2 = fc 2 such 
that 

k 2 = -k, (117) 

In terms of the (axial) dispersion relation uj (k), the Eq. (|117|l boils down to a simple 
definition of axial spectral symmetry 

u [k x , ky,k) = Co {k x , ky, fc) , 

where ki = k and fc 2 = — k constitute a pair of reciprocal real wave numbers related 
to given w and k x ,k y . In the case <|107[) of four propagating eigenmodes, there will be 
an additional pair k^ and ki = — k% of reciprocal wave numbers. 

Now assume that while k\ and fc 2 = —kx are real, the remaining wave numbers 
fea and ki from the set H116(l are complex, which constitutes the case H109JI of two 
propagating and two evanescent eigenmodes. In such a case, in addition to Eq. (|117(l . 
the relation 1(116(1 together with ((106(1 yields 

k 4 = -k 3 = k* 3 , (118) 

or cquivalently 

Re ki — Re k 3 = 0, tt/L; Imfc4 = — Imfc3. (H9) 

In Eq. I|l 19(1 we took into account that k = k + 2ir/L. 

Consider now the case 1(111( 1 of a frequency gap, where all four eigenmodes are 
evanescent. The relations ((11611 and 1(1 U6|) allow for two different possibilities. The 
first one is similar to that of Eq. 1 (118(1 



k 2 = —kx = k*, 



or, equivalcntly, 



(120) 

Re fc 2 = Re k% = 0, tt/L; Imfc 2 = — Imfci, (121) 
Re ki = Re k% — 0, tt/L; Imfc4 = — Imfc3. 



In the above situation, the four complex wave numbers split into two reciprocal pairs 
kx, fc 2 and &3, ki of the conjugate values. The other possibility is 

kx = -fc 2 = feg = -ki, (122) 

or, equivalently, 

Re kx — Re k% = — Re ki — — Re ki, Im kx = — Im /C3 = Im fc 2 = — Im ^4. (123) 



5.4-1- Spectral asymmetry If the sufficient condition 1(113(1 for the axial spectral 
symmetry is not in place, then we can have for given u> and k x , k y 

{kx, k 2 ,k 3 ,k 4 } {-kx,-k 2 ,-k s , -k 4 }, (124) 

which implies axial spectral asymmetry 1(67(1 . The relation ((106(1 . being a direct 
consequence of J-unitarity ((91(1 of the transfer matrix Tl, remains valid. 
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5. 5. Electromagnetic energy flux in stratified media 

5.5.1. The J -scalar product For future reference, consider the following scalar 
product involving the J-matrix i|86|) 

(*i> ^2) = El x H 2y - E\ y H 2x + Hl y E 2x - H*[ x E 2yi 

which will be referred to as the J-scalar product. Given the importance of the above 
quantity, hereinafter, we will use the following special notation for it 

[*i,* 2 ] = (*i,J* 2 ). (125) 

The J-scalar product (|125[1 is invariant under the following transformation involving 
an arbitrary J-unitary matrix T 

[T*i, T* 2 ] = *2] for any *iand * 2 . (126) 

The relation (|126fl can also be viewed as a criterion of J-unitarity of a matrix T. This 
relation is similar to that involving the regular scalar product (^1, ^ 2 ) and a unitary 
matrix U 

{U^i, U^ 2 ) = (*i, * 2 ) for any *iand ^ 2 . 

Let (z) and (z) be two arbitrary solutions of the time-harmonic Maxwell 
equation l(%Tjl. The equality l|126|l together with the definition (f%%|) of the transfer 
matrix yields 

[* {z\ , 9 j (z)] = [T (z, 0) * (0), , T (z, 0) 9j (0)] = [* (0), , (0)] , (127) 

which implies that the J-scalar product [^f (z) i , ^ j (z)] does not depend on the 
coordinate z. 

Consider now the J-scalar product 

(128) 

of two eigenvectors $i and <!>., of the transfer matrix Tl 

T L $ i = ( i $ i , t= 1,2,3,4. (129) 
The J-unitarity (|126[) of T L implies that 

[T^Ttfj] = CiCj *j] = *j] > (130) 
which, in turn, yields the following important relation 

[4^1=0, if CC^i, (131) 

or equivalently, 

[$i,<&i] = 0, if k^k*. (132) 

In particular 

[$i,$i]^Q, only if & = (133) 
which means that <£>i in H133(l should be a propagating Bloch mode. 
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5.5.2. Energy flux in stratified media The real- valued energy flux (the Poynting 
vector) associated with a time-harmonic electromagnetic field is 

S (r) = [ReE(f) x ReH (r)] = i Re [E* (r) x H (r)] . (134) 

Substitution of the "tangential" Bloch representation (|80J) for E (r) and H (r) in Eq. 
ifTMjl yields 

S(r) = S(z) = ^Re E* (z) x H (z) , (135) 

at fixed u> and k x ,k v . Eq. I|135|) implies that in a stratified medium, at fixed u and 
k x ,k y , all three Cartesian components of the energy flux S (r) are independent of the 
tangential coordinates x and y. A simple energy conservation argument shows that 
the normal component S z of the energy flux does not depend on the coordinate z 
either, while the tangential components S x and S y may depend on z. Indeed, in a 
lossless stratified medium we have, with consideration for Eq. I|135(l 

v • s (f) = d z s z (z) = o, 

which yields that at fixed u> and k x , k y 

S z (r) =S Z = const. (136) 

By contrast, the tangential components of the steady-state energy flux are dependent 
on the z coordinate 

S x (r) = S x (z) , S y (r) - S y (z) . (137) 

Hereinafter, the normal component of the energy flux will be referred to simply 
as the energy flux, unless otherwise specifically stated. It also will be denoted as S, 
rather than S z . 

The explicit expression for the normal component of the energy flux l|135|) can be 

recast as 

s = \ (KH y - e;h x + e x h; - e v h*) = ~ [9, 9] , (138) 

where the J-scalar product = (*, J*) is defined in Eq. (|T25)l . The fact that S 

in Eq. (|138|l is independent of z implies that 

S = i (z) , * (*)] = i *] , where $ = * (0) . (139) 

Eq. i|139fl can also be viewed as a direct consequence of J-unitarity l|91|l of the 
transfer matrix. Indeed, from the definition (|88(l of the transfer matrix we have 

V(z) =T(z,0)*(0) =T(z,0)$. (140) 

Substituting ljTIU|) into (|HS|> yields 

S = i[* («),*(«)] = i[T(z,0)$,T(z,0)$]. 
Taking into account the property (|126|l of a J-unitary matrix, we again arrive at Eq. 
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5.5.3. Energy flux in periodic stratified media The direct relation l|139f) between 
the J-scalar product [iff, iff] and the energy flux S at fixed u> and k x , k y allows us to 
make some strong statements regarding electromagnetic energy flux in periodic layered 
media. 

Let us start with the simplest case of a single Bloch eigenmode. Eq. 1(139(1 
together with 1)133(1 shows that only a propagating mode can transfer electromagnetic 
energy 

Si = ^ ^ only if hi = k*. (141) 

A single evanescent eigenmode always has zero energy flux 

Si = [$,, = 0, if fc^fc*. (142) 
Let us turn to the case of a superposition 



i=l 

of different Bloch eigenmodes with fixed w and k x ,k y . In such a case, the energy flux 
is 

1 1 4 

S= £ 2 .? <|, - ,|> - ( 143 ) 

Taking into account Eqs. ((132(1 we can draw the following conclusions: 

1) The contribution Si of each propagating eigenmode to the total energy flux is 
independent of the presence or absence of other Bloch eigenmodes with the same 
uj and k x , k y 

i=l i=l 

where the summation runs over all propagating eigenmodes. The number of 
propagating modes can be 4, 2, or 0, depending on which of the cases 1(111(1 . 
1(109(1 . or 1(111(1 we are dealing with. 

2) The contribution of evanescent Bloch eigenmodes to the energy flux depends on 
their number. 

(a) In the case ((1U9(I of two evanescent modes $3 and $4 we have 

S = Re (0304 [$3, $4]) , where fc 4 = fcg, (145) 

which implies that only a pair of evanescent modes with conjugate wave 
numbers can produce energy flux. The respective contribution ((145(1 is 
independent of the presence of propagating modes $1 and $2- In accordance 
with Eq. ((142(1 . a single evanescent mode, either <I>3 or $4, does not produce 
energy flux on its own. 

(b) In the case ((111(1 of four evanescent modes we have 

S = Re (a*a 2 [$1, $2]) + Re (0304 [$3, $4]), where k 2 = k\ , fe 4 = fcg , (146) 

which implies that either of the two pairs of evanescent modes with conjugate 
wave numbers contribute to the energy flux independently of each other. 
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5.6. Scattering problem for a periodic semi-infinite stack 

In this final subsection we outline the standard procedure we use for solving the 
scattering problem of a plane monochromatic wave incident on the surface of a periodic 
semi-infinite stack. 

In vacuum (to the left of the semi-infinite slab) the electromagnetic field f y ( z ) 
is a superposition of the incident and reflected waves 

*y (z) = i& x (z) + V R (z), z< 0, (147) 

where the indices I and R relate to the incident and reflected beams, respectively. At 
the slab boundary we have 

$ v (0) = *j(0) + * fl (0). (148) 

The transmitted wave (z) inside the periodic semi-infinite slab is a 
superposition of two forward Bloch eigenmodes 

^ T (z) = *i (z) + $2 (z) , z> 0. (149) 

The eigenmodes (z) and \&2 (z) can be both propagating (with u > 0), one 
propagating and one evanescent (with u > and Imfc > 0, respectively), or both 
evanescent (with Imfc > 0), depending on which of the three cases l|39(l . 1|41|1 . or (|40ll 
we are dealing with. 

Assume now that for a given frequency ui, the Bloch eigenmodes are found, which 
can be readily done in the case of a periodic layered array. Using the standard 
electromagnetic boundary conditions 

*r(0)=*/(0)+*ji(0), (150) 

one can express the reflected wave "f^ and the eigenmode composition of the 
transmitted wave ^x, in terms of the amplitude and polarization of the incident 
wave ^j. This automatically gives the electromagnetic field distribution (z) inside 
the slab, as a function of the incident wave frequency, polarization, and direction of 
incidence. 

The transmittance and reflectance coefficients of a lossless semi-infinite slab are 
defined by the following expressions 




(151) 



where (jSij , \&Rj an( i \^ T ) arc norma l components of the energy flux of 
the incident, reflected, and transmitted waves, respectively. Knowing the value of the 
transmitted wave ^S?t or reflected wave VE^ at the slab boundary, one can immediately 
find the respective energy flux and, thereby, the transmittance/reflectance coefficients 

cm. 

The above-outlined standard procedure was used in all our numerical simulations. 
It applies both to the case of normal and oblique incidence. In the latter case, the 
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explicit expressions for the column vectors fyj and ^>r in I|147I150[) are 









Er, x 




E I,V 

Hi,x 




E R,y 
Hr,x 



(152) 



where the complex vectors Ex, Hi and Er,Hr are related to the actual 
electromagnetic field components E/,H/ and E.r, Hr as 

Ej = e ^(»^+«^) j g / ( z ) t H 7 = e^^+^Hi, (153) 
= e ^ {n * x+n ^E R (z) , H R = e ^i n ^+ n yV)H Ri (154) 

as prescribed by Eq. (|80|l . Here n is the unit vector in the direction of light propagation 

/ c c c \ 

for incident beam: n = fti = ( — k x , —k v , —k z ) , (155) 

\LU UJ Ul ) 



/ c c c \ 
for reflected beam: n = Hr = ( — k x , —k v , k z ) , (156) 

\UJ UJ to J 



where 



c c 
— n 2 = — 

UJ UJ 



Note that the tangential components of the unit vector n of the incident wave are 
the same as those of the reflected wave. The electric and magnetic fields of a plane 
monochromatic wave in a vacuum are uniquely related to each other 

H = n x E. 

The same relation holds for the complex vectors Ei, Hi and Er, Hr defined in l|153fl 
and ^ 154(1 . namely 

Hi =ni x Ei, H R = n R x E R . (157) 



6. Matrix of reflection coefficients of a semi-infinite periodic stack 

This and the following sections are devoted to a rigorous mathematical analysis of the 
scattering problem for a plane monochromatic wave incident on a periodic semi-infinite 
stack. We focus on vicinities of stationary points of the electromagnetic dispersion 
relation and our goal is to develop an asymptotic analytical description of the frozen 
mode regime. Not only that would allow to rigorously prove the physical results 
presented earlier in this paper, it would also provide a better understanding of the 
very essence of the frozen mode regime. The major part of the following analysis 
is a perturbation theory of degenerate non-diagonalizable matrices. Specifically, we 
refer to the transfer matrix Tl, which develops a nontrivial Jordan block at any 
stationary point of the dispersion relations. The latter circumstance implies the 
existence of diverging non-Bloch eigenmodes, which usually do not contribute to the 
transmitted wave inside the semi-infinite photonic slab and, therefore, do not 
affect the scattering problem at hand. Yet, there are two important exceptions. The 
first one is the stationary inflection point (|15|) . where not only the linearly diverging 
Floquet eigenmode dominates the transmitted wave, but it also produces a finite 
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energy flux inside the periodic medium. Another exception is the degenerate band 
edge where the respective linearly divergent non-Bloch eigenmode, although 

dominant, does not contribute to the energy flux and, therefore, does not effectively 
transform the incident radiation into the slow mode. 

The rest of the paper is organized as follows. In this section we re- formulate 
the scattering problem for a lossless periodic semi-infinite stack, introducing basic 
notations and definitions. In the following sections we develop a perturbation theory 
for degenerate non-diagonalizable 4x4 matrices and apply this theory to the transfer 
matrix Tl and, thereby, to the scattering problem. Special attention is given to the 
comparative analysis of different stationary points of the electromagnetic dispersion 
relation, such as a photonic band edge, a stationary inflection point, and a degenerate 
band edge. To simplify the rather cumbersome mathematical expressions of the 
following sections, we will use the following new notations for the quantities already 
defined earlier 

xi — > x, x 2 -> y,x 3 -> z, 

A — > -JM, T^T L , 
c 

k = (ki,k 2 , ks) — > (ck x ,ck y , ck z ) . 
Observe that the 4x4 matrix 

A = -JM (158) 

c 

is Hermitian, while the related Maxwell operator M defined in Eqs. (|81I86|) is J- 
Hermitian. 



6.1. Basic definitions 



A periodic semi-infinite stack is defined in terms of the related matrix function A (2:3) 
satisfying 

A (x 3 ) = Const, -00 < x 3 < 0; A (x 3 + L) = A (x 3 ) , < x 3 < 00. (159) 
In vacuum, the Hermitian matrix A (x 3 ) defined in (|158f) and has the form 



A(x 3 ) = AM = 



a(°) 
a<°) 



, a<°> = - 



kik 2 



kik 2 
uo 2 -k\ 



(160) 



The above expressions immediately follow from Eqs. H83[) and 1)158(1 . The tangential 
component k r of the the wave vector k is related to its normal component k 3 by 



V^ 2 -k2 ^^ 2 -(fc2+fc 2 2 ). 



Let us introduce 



and notice that 



.12 



-1 

1 



J = 



-h 

h 



: 32® 32, 



32* 



(0) 



-J2a 




(0) 



-1 

1 



32a 



(0) 



: J2 ® 



J2a 



(0) 



(161) 



(162) 



(163) 
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Recall the basic properties of the tensor product operation: if A and B are square 
matrices and u and v are vectors of related dimensions then 



[A <g> B](u®v) = {Au) ® (Bu) , (ui ® «i,U2 <8> U2) = (fij "2) («i,«2) 



(164) 



Suppose now that we know the set of eigenvectors and eigenvalues for two square 
matrices A and B, namely 



Auj = XjUj, Bv m = n m v m . 

Then ((Pol}) and imply 

[A (g) B] (uj ®!) m ) = \jfj, m Uj <g> u„ 



(165) 



(166) 



Using (|166|) and the tensor product representation (|163|l for JA^ we can find its 
eigenvectors and eigenvalues as follows. First, we find that 



j 2 u± = ±iu±, U± = 
v± = v± (lo, k r ) = 
Notice also that 



1 
1 



±i 
1 



; j 2 a. (0) v ± = ±ik 3 v ± , 



(167) 



±io;A;3 + fciA;2 

1 



1 2 (kl + kj) k 3 uj 



(v+J 2 v-) = (v-,j 2 v+) = 0, (v±,j 2 v±) = ±i 
(v±,v±) = 

1 



(168) 



(fc| + k\ 



2 + (cofc 3 ) 2 + (kik 2 ) 2 



(v T ,v ± )- 2 



^3 + ^2 



2 (fcg + fc|) fc 3 W 
2 

+ 1 



±iwfc3 + A:ifc 2 



(u±,u±) = l, (u T ,u±) = 0; (u T ,j 2 u±)=0, (u±,j 2 u±) = ±i. (169) 

Using the tensor product representation i|163|) for JA^ ^ and (|166fl and (|167|l . I|168|l 

we obtain 



JAW 



Zf = ±k 3 Zf 



JAM 



Z±=±k 3 Z±, Zf = Zf(u,K). 



where 



Zf = u- (g) v + , Zf = m+ <E> U- , #1 = u+ ® v+ , Z 2 = U- tg) V- 



The component representations for Z ■ are as follows 



i(fcg+fcl) 

ifc 3 w+fei fe 2 

— i 



ifc3^ + fel &2 
1 



' Z 2 



i(fcl+fcl) 

— ifc3CJ+/cifc 2 

i 

— ife3CJ + /cifc2 

1 



(170) 



(171) 



(172) 
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ik^u+ki k 2 

i 

fcg + ^2 

ife3^+fei k 2 
1 



V2 7 , 



ifol+fcjl) 

— i 

— ifc30!+&l fc2 
1 



(173) 



(174) 



(175) 



(176) 



Observe that (|ToI|> . ((T7T|l . imply 

[Z+, Z+] = (Z+, JZ+) = («_ <g> «+, [j 2 ® ia] «+ ® v-) 
= (it_ ® u + , j 2 u + ® J2«-) = i («- ® <g> 

= (u+,y'2«-) = o, 

[ z i> z t] = ® ®j2f+) = -i(w-,u-) (u+, j 2 u+) = 1, 

\ z \i z x \ = ® <g> j 2 u + ) = -i(u_,u+) (i>+, j 2 u + ) = 0. 

Carrying out more evaluations similar to 1174|) . i|175|) we get 

[Zf, Z%\ = 0, [Z±, Z±] = ±S jm , j, m = 1,2, 

where (5j m is Kronecker symbol. The relations (|176l) show that the system of 4 vectors 
Zj , j = 1,2 is flux-orthonormal in the sense that it is orthonormal with respect to 
the flux form * 2 ] = J^i)- 

Consider now the scalar products of 4 vectors Z z > , j = 1, 2: 

(Z+,Z+) = («_ <g> ®t; + ) = (u_,u_) («+,«+) = U ^0 V , (177) 

{ z i: z i) = (u + ®v + ,u + ®v + ) = (u + ,u + )(v+,v+) =/3 U)kT , 

(Zf,Z^) = {u^<E)v + ,u + <E)v+) = (u^,u+){v+,v + ) = 0, (178) 
(Z^ , = (u_ ® v+, u+ <8> w_) = (u+, u_) = 0. 

Carrying out evaluations similar to (|177f) . I|178|) we get the following complete set of 
equalities: 



{Z?,Z±) = Uikr 8 jm , (Z?,Z%) = 0, j, m = 1, 2. 



(179) 



(fc 3 2 + fc 2 2 ) + M&) + 0W 



showing that the system , j 
orthonormal. 



2(fc 3 2 + fc 2 2 )fc 3 

1,2 is orthogonal though, evidently, it is not 



The set of equalities l|176f) and (|179f) show the system of vectors Z. , j = 1,2 has 
a property that both the forms, namely, the EM density form (the scalar product) and 
the flux form, become diagonal if it is chosen to be a basis of the space C 4 . Another 
advantage of choosing Z^, j = 1, 2 to be a basis is that in this basis the flux balance 
equality for relevant modes takes its simplest form as in the classical scattering theory 
(see (EHJ, (|215|) . In fact, the latter is our primary motivation. 

For the periodic semi- infinite stack with A (x 3 ; u>) we have the following equation 
defining its eigenmodes \& (£3) at the frequency u 



(x 3 ) = iJA (x 3 ) * (x 3 ) , -00 < x 3 < 00. 



(180) 
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The eigenmodes of the periodic semi-infinite stack are the ones corresponding 
to an incident wave which propagates from — oo to oo, then it is partially reflected 
by the interface at x 3 = and partially transmitted into the dielectric substance in 
< X3 < 00. We refer to such eigenmodes as relevant eigenmodes and denote the set 
of all relevant eigenmodes by S T = S T (uj) = S T (u, k T ). 

The two extended eigenmodes \?i (0:3) and \t/ 2 (£3) describing the standard 
scattering problem satisfy the following relations in the air, 2:3 < 0, 



vl/i (as 3 ) = e^Zj 
^2 (x 3 ) 



e 



ik 3 x 3 7 + 



— lk 3 X 3 

-ik 3 x 3 



[Pll Z l +P21 Z 2 ] > 



(181) 



where the matrix of reflection coefficients 

P = Pu>,U.r z 



p n (u),k T ) p 12 (u),k T ) 
p 21 (uj,k T ) p 22 (uj,k T ) 



(182) 



carries the information about reflection properties of the slab. Its entries can be called 
reflection coefficients. 

The set S T of all relevant eigenmodes \& (x 3 ) happens to be a two-dimensional 
linear space. For every fixed real a it is uniquely determined by the two-dimensional 
space S T (a, to) = S T (a; ui, k T ) of the values VP (a) as \t runs over S T (k T ), i.e. 



S T (a; w, k T ) = {* (a) : * £ S T (k r )} . 



(183) 



More precisely, all possible relevant eigenmodes are described by solutions to the 
following Cauchy problem 



^3* (£3) = iJA (13) * (x 3 ) , *(a) = $e5 T (a, to) , -00 < x 3 < 00. 



(184) 



The two-dimensional space <S T (a, lu, k T ) provides a convenient way to describe and 
parametrize the relevant modes. For instance, assuming that we know <S T (a,u, k T ) 
let us pick any $ £ <S T (to, k T ) and find values of the eigenmode vp (a; 3 ) in the air. 
The eigenmode >J/ (aj 3 ) can be represented as the following linear combination for 

—00 < X3 < 00 : 



*(z 3 ) 



[a+Z+ + a+Z+] 



* (0) = Z+ + a+Z+ + o^Zf 
where evidently the two pairs of coefficients 



[ajZf+aaZa-] 
Z 2 -<Se5 T (0), 



(185) 



and a = 



a. 



(186) 



are respectively related to the incident and the reflected waves. As is commonly done, 
we choose arbitrarily the incident wave by picking the vector a + and then finding the 
the reflected wave as the vector a~ using the relations H181|) . (|182fl and (|185[) by the 
following formula 



a 



pa + , a 



P11 P12 
P21 P22 



(187) 
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Observe then that the matrix of reflection coefficients p can be viewed as the following 
mapping relating the incident wave a + to the reflected wave a~ 

p:a + ^a~. (188) 

Notice also that the reflection and the transmission coefficients r (a + ) and t(a + ) 
corresponding to the incident wave a + are defined by the formulae 

r> K) = = ^ t2 ^ = 1 " ^ (189) 

|a+($)| \a+\ 

It follows from 1|> that the space 6>t (0;cJ,k T ) has the following representation in 
terms of the vectors Zj and the reflection coefficients Pj m : 

S T (0; w, k T ) = Span {(Z+ + Pll Z{ + p 21 Z 2 ) , (Z+ + p 12 Z^ + p 22 Z 2 )} . (190) 

The relation (|190[) shows that the space Sj> (0;w, k T ) is uniquely determined by the 
matrix p^ kT . We show in the following subsection that the matrix p^ k ^ is uniquely 
determined and can be constructed based on the space S T (0;w, k T ). 

6.2. Basic properties of the space of relevant eigenmodes 

Let us consider now basic properties of the two-dimensional space (0; ui) = 
<St(0;u, k T ) suppressing in the notation its dependence on k r . Notice first, that 
the space <St (0; uj) has the fundamental property that it is always nonnegative with 
respect to flux form in the sense that 

$] > for any $e5 T (a; u)) . (191) 

The property l|191|) indicates that the modes related to S T (a; uj) must transport energy 
in the chosen direction. 

It is a well known result of spectral theory that no eigenmode '3/ (£3) can grow 
at infinity faster then polynomially. In particular, an eigenmode can not grow 
exponentially as 13 — > 00. Since \& (x^) is a solution to 118Ufl it must be a linear 
combination of eigenmodes of the infinite periodic stack with the relevant periodic 
A(x3;lo) on the interval < £3 < 00. Consequently, such a linear combination can 
not include evanescent modes growing exponentially as S3 — > 00. Additionally, the 
above-mentioned linear combination cannot include backward propagating eigenmodes 
(those with negative group velocity). Notice that, the related properties of S T (0;uj) 
can be characterized by the spectral properties of the transfer matrix T (oj) . For 
instance, in the case when all eigenmodes are propagating and have different wave 
numbers as described by H108(l . the space S T (0; u>) is the span of those two eigenvectors 
${ 1 = (lj) and <I>£ 2 = (uj) that have positive fluxes, i.e. 

5 T (0;o;)=Span{%,$ e2 }, where [$f 2 , * Ca ] > 0. (192) 

Hence, there are exactly two eigenvectors having positive fluxes. 

In the case l|110|) when there are two propagating and two evanescent modes, 
6> T (0; ui) is the span of the two eigenvectors (uj), having a positive flux, and $^ (ui) 
with |(| < 1, i.e. 



5 T (0;w)=Span{^ 1 ,$ c }, where \^\ = 1, , > and |£| < 1. (193) 
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Finally, in the case (|112f) when all modes are evanescent, S T (0; uS) is the span of the 
two eigenvectors (ui) and $<^ 2 (u>), i.e. 

5 T (0;w) = Span{$ fl ,$ C2 }, where |Cil,|C 2 l<l- (194) 

If for a certain frequency loq the transfer matrix T (wo) has a non-trivial Jordan 
block, then the space <S T (0-ujq) can be defined as the following limit 

S T (0;w ) = lim S T (0;u), (195) 

UJ — >LUo 

where it is assumed that for lo ^ luq the matrix T (w) is diagonalizable and (0; w) 
is a well defined two-dimensional space. The limit l|195|) uses a distance d between a 
two subspaces Si and 5*2 defined by the formula, 08] j Section IV, §2, 

d(Sx,S 2 ) = max {S (Si, S 2 ),S(S 2 , Si)}, S(S 1 ,S 2 ) = sup sup ||u - v|| . (196) 

ueSi:||u||=l v es 2 

Hence, the limit relation in (|195fl is interpreted as 

lim d(S T (0;ui),S T (0;ui )) = 0. (197) 

u? — »aJo 

The distance d(Si,S 2 ) defined by l|196|> measures the "aperture" or "gap" between 
the subspaces Si and S 2 . It has the following important property, 08], Section IV, 
§2, Corollary 2.6, 

d(Si,S 2 ) < 1 implies dimSi = dim£ 2 . (198) 

The property (|198|l implies that if the limit (|195|) exists then the dimension of the 
space S T (0;wo) must be 2 since dim<5> T (0;w) = 2 for u ^ u>o- 

One can also verify that the limit relations l|195|) , (|197|l can be conveniently recast 
as a limit relation between orthogonal projections onto spaces <S T (0;u>). Namely, if 
we introduce 

Ps to be the orthogonal projector on the space S, (199) 
\\P S \\ = sup ||Ps$|| , where ||$|| is length (norm) of ||$|| . (200) 

then (|195|l . (|197|l are equivalent to 

lim \\Ps T (^) - p s T (0:u )\\ =0. (201) 

UJ — >U)Q 

Notice, that the relation (|201|l is equivalent, in turn, to the relation 

lim Ps T (o ; a,)$ = Ps T (o-, Uo )$ for any $ e C 4 , (202) 

UJ — *UJq 

and the relation (|202|l is equivalent to 

lim P St (0;u,)$ = Ps T (0;^ )$ = $ for any $e5 T (0;w q ) . (203) 
Notice now that for every vector $ G S T (0; u;o) we can define a family of vectors 

$H = Ps T (0;^)*eS T (0; W ) (204) 
converging, in view of (|203|l . as u> — > Wo to the vector $ (wo), i.e. 

lim $ (u>) = $ (o>o) = for any $e5 T (0; Wq) . (205) 

UJ — >u>o 
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In this section we look at the basic properties of the matrix of reflection coefficients 
P = Puj k T as defined in l|18ip . I(182|) . (|187|l . and its relation to the flux quadratic form 
[■, ■], and the space S T (0;u>, k r ). 

Observe that inserting X3 = into (|185|) yields 



f (0) = $ = $+ + $± = of ($) + af ($) ^2 



(206) 



The equality 1(206(1 indicates that the numbers ($) are the components of the vector 
$e5 T (0; w) with respect to the basis 

{Z+,Z+,Zi,Z^}, (207) 

and they are determined by the following formulae. Let us introduce the two- 
dimensional subspaces of C 4 : 



Z+ = Span {Z+, Z+) , Z~ = Span {Zf , Z^ } , 
and the respectively orthogonal projections: 

7r + and 7r~ are respectively the orthogonal projections on Z + and Z~ 
In view of 1179(1 we have the following representations for ir±: 



7^$ = 



1 



[(zf 1 <s>)zf + (z±,<s>)z±], $e 



and, hence, 



$ = + $ ± = 71-*$ = ai ($) + ($) 
1 



/3a 



\ziM 
, 2 



(208) 
(209) 

(210) 

(211) 
(212) 



U?,k r 



(fc^ + fc 2 y + (a;fc 3 r + (fcifc 2 r 
2(fcf + fc 2 )fc 3 



&3 = 



k?. 



Observe, in particular, that for <i> G i5 T (0;a;, k T ) i/ie equalities \21ty) provide relations 
between the value 4> of the mode at X3 — and the coefficients a (<&) for the relevant 
incident and reflected waves. 

Another simple fundamental fact is that the two-dimensional vector a + ($) can 
take any prescribed value from C 2 , i.e.: 



{«+($) : $e S* T (0;w,k r )} = C 2 , a* (*) 



(zt,*) 



(213) 



The relation (PTSfl can be considered as another fundamental property of the space 
S T (0;uj,k T ). 

Using the coefficients a ($) and 1(176(1 we get the following representation for 
the flux of the mode described by $ 



[$,$] = [$+,$+] - [$",$"] = \a + ($)| 2 - I a - ($)| 2 , 
$ = afZ^ + a\Z^ + a\Z^ + a^Z^ e S T (0;cj,k T ) . 



(214) 
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The above equality reflects the fundamental energy flux balance of the classical 
scattering theory in its the simplest form: 

\ a+ (^ ) )| (Incident wave flux) - | a - ($)| (Reflected wave flux) (215) 
= |a + ( < 1 > )| 2 — \a~ ($)| 2 (Transmitted wave flux). 

The fundamental property (|191f) of the non-negativity of the flux on <S T (0; u, k T ) 

can be recast as 

| a~ ($)| 2 = [$~,$~]< \a+ ($)| 2 = [$+,$+] for any$GS T (0;w,k T ), (216) 

indicating the physically transparent fact that the flux of the reflected wave cannot 
exceed the flux of the incident wave. 

Combining now the relations (|187fl and (|216[) and recalling that a~ ($) = 
pa + ($), gives 

\pa + ($)| 2 < \a+ ($)| 2 for any a + ($) (any $ <E 5 T (0; k T ) ), (217) 

which, in turn, together with 1213|) implies 

pt p < J 2 . (218) 

The matrix inequality (|218|l signifies the fact that for any or any incident wave a + , 
the reflection coefficient r (a + ) does not exceed 1, i.e. 



r2 (a+) = \^mf = \p« + (*)i 2 = k4 = K> ^ pJ f , < L (219) 

V 7 |a+($)| 2 |a+($)| 2 |a+| 2 (a+,a+) 
As to the EM energy density using (|179|l we get 

($,$) = /?, (|a+| 2 + |a-| 2 ), (220) 



(kl+klY + {Ljk 3 f +{k X k 2 ) 



2 



^> k = 2lfcf+fc|)fc^ ' fc3 = 

$ = + + aT/Zf + Z% G 5 T (0; w, k T ) . 

Having the energy flux balance in the form (|214|l was our primary motivation 
for choosing the vectors Zj, j = 1,2 as a basis in C 4 . We also want to remind the 
reader that the vectors Z^, j = 1,2 reduce both the EM energy density (the scalar 
product) and the flux quadratic forms to their diagonal form, as follows from the set 
of equalities and HUB- 

Let us look now at the limit case p> p = I2 for which, according to l|214|) . 

[$,$] = \a + ($)| 2 - | a - ($)| 2 = for all a + G C 2 , a~ = pa + , (221) 
$ = aXZ^ + a\Z^ + af-Zf + a% Z% G S T (0; w,k r ) . 

If we denote by ©o («/) the set of spaces on which the flux is identically zero, then the 
relations (|221[l and l|190|) imply that 

S T (0;w,k r ) = (222) 

= Span{(Z^ +p n Zi +P21Z2) > ( z 2 + Pu z i + P22 z 2~)} e ®o (J) 
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or, in other words, all vectors of the space S T (0;w, k T ) have zero flux. On the other 
hand, if S T (0;w, k r ) £ <3 (J) then (|221[1 holds implying p 1 p = I 2 . Consequently, the 
property that the slab has complete reflection p^ p = I 2 is equivalent to the property of 
having zero flux for all relevant modes, or, symbolically, 

p*p = h is equivalent to [$, $] = for every $e5 T (0; u, k r ) , (223) 



or, in other words, 

the reflection coefficient r (a H 



1 for every a + £ C 2 is equivalent 



(224) 



to [$ (a + ) , $ (a + )] = for every a + £ 



Therefore, to establish the state of complete reflectance it is sufficient to verify that 
the fluxes of all relevant modes are zero. 

Observe also that as a consequence of {2236 . {22$ we have 



if there exists $e5 T (0; u, k T ) such that [$, $] ^ 0, then p ] p ^ I 2 , (225) 

or, in other words, 

if there exist $e5 T (0; to, k T ) such that [$, $] ^ (226) 
then for almost all a + £ C 2 : the reflection coefficient r (a + ) < 1. 

To establish a representation for the matrix p in terms of the space S T (0; u>, k T ) 
let us pick any two linearly independent vectors $1 and $2 in 5 T (0;w, k T ). Then, 
since 5 T (0;w, k T ) is a two-dimensional space, we have 



S T (0;w,k T ) = Span {$1, $2} ■ 



(227) 



Having the basis {$i,$ 2 } of <S T (0;w,k T ) we introduce the related component 
representation 



4> 



$ = (p^i + ip 2 $ 2 £ S T (0; ui, k T ) 



and 



7^$ = a* ($) = Q ± $, 
Q± = [Z+Z+] t [$ 1 $ 2 ] 



(228) 



(229) 



i 









{zfM) . 



Observe now that the relation a = pa + together with l|229|) implies 
where 

$ = $ (a + ) = [Q + ] _1 a+, $ (a+) = [$i$ 2 ] $ (a + ) = [$1*2] [<3 + ] " + 
The relations l|230|) yields, in turn, the following representation for the matrix p 



Q- [Q H 



(Z+,*!) (Z+,$ 2 



(230) 
(231) 

(232) 
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Notice that the inequality (|216[) together with i|213|) implies 

[Q-^Q- < [0 + ] f + , (233) 

which is an alternative form of the inequalities l|216|) . (|218|) and l|219|) . Using (|219|) . 
I|229|) and i|232|) we get the following representation for the reflection coefficient 



{q-[Q+} 1 a+,Q-[Q+] \ 



r 4 K) = j-T]3 S (234) 

|a+| 

^KW)^^^^. (235) 
V 1 " |a+(*)| (Q+*,Q+$) 

Observe also that 1)22 l|l and 12341) yield the following expression for the flux associated 
with the incident wave described by a + 

[<D (a+) , $ (a+)] = (1 - r 2 (a+)) |a+| 2 (236) 



l«-| 2 



The formula l|236|) can be recast as the following representation for the transmission 
coefficient t = t (a + ) defined by l|189l) 



e (a+) = 1 - r 2 (a+) = [$ ^ ' $ 2 ( " +)] , $ (a+) = [«x* a ] [Q+] «+. 



7. Transfer matrix at and near a point of degeneracy 



(237) 



Let us recall first the definition of degenerate points including inflection ones. A n - 
degenerate point fco of a dispersion relation u> (fc) is defined as a point at which the 
following relations holds 

d k u (fc ) = d 2 k cu (k ) = ■■■ = d r ^uj (k ) = 0, dluj (k ) ? 0. (238) 

In particular, an inflection point fco is a 3-fold degenerate point if 

J (ko) = oj" (k ) = 0, J" {k ) ^ 0. (239) 

Hence, if fco is a n-degenerate point we have 

u (fc) = uj (fc ) + ^ [ fco) (fc - fc )" + O ((fc - fc )" +1 ) , fc -» fc . (240) 
In particular, if fc is an inflection point then 

w (k) = w (fc ) + ^1^1 (fc - fc ) 3 + o ((fc - fco) 4 ) , fc -> fco- (241) 
To study the behavior of the transfer matrix Tl near loq we introduce 

T (v) — T (ujq + v) , f = w - oj . (242) 



Slow light in photonic crystals 52 

We assume the dependence of T (v) on v to be analytic in some vicinity of v = 0. 
In our further analysis we use well known statements from the analytic perturbation 
theory for matrices and their spectra, |48| . 

To find the spectrum of T (y) we consider the characteristic polynomial ^t{v) (C) 
and the related characteristic equation 

At» (0 = det (T {v) - Qh) = 0, C = e lk , (243) 



where I4 is the 4x4 identity matrix, and k is the quasimomentum. Eq. I|243[l is the 
dispersion relation, namely it relates to every frequency v — to — uj four values of C 
or, equivalently, four values of the wave number (see Eq. H1U4|) '). 

Since T {v) is 4 x 4 matrix the equation (|243[1 can be written as 

Ar ( ,) (C) = C 4 + b 3 (1/) C 3 + b 2 (1/) C 2 + h (i/) ( + b (y) = 0, (244) 

where the complex valued functions bj (v), j — 0, 1, 2, 3 are analytic in v in a vicinity 
of v = 0. 

For the frozen mode regime to occur at the frequency v — 0, i.e., lu = ujq, the 
spectral decomposition of the transfer matrix T (0) must have a Jordan block of rank 
n > 2 with an algebraic eigenvalue Co- I n this situation the characteristic polynomial 
At-(q) (C) takes the following special form 

A r( o) (0 = (C - Co)" Qn (0 , (245) 

where 

Qn (C) — C 4 ~" + - ■ ■ is a polynomial of the degree 4 — n such that Q n (C ) 7^ 0. (246) 

It is an additional property of the transfer matrix T (y) that 

ICol = 1, (247) 

where the eigenvalues Co = e%ka 1S n-degenerate and it corresponds to a Floquet mode. 
Because of this degeneracy at v = 0, the perturbation theory, 0H]i Section II, classifies 
the point v = as an exceptional one, and the dependence Co iy) is described by the 
Puiseux series of the form 

Co M = Co (l + uiv 1/n + <*2V 2/n + •••)• (248) 

The corresponding eigenprojectors can be singular. In fact, in our case they are 
singular. 

If the characteristic equation (|243|l takes the special form l|245f) near v = then 
T {y) can be reduced and represented as follows 



T{v) = Q{v) 



7» 
W{v) 



Q- 1 {v) , (249) 



where Q (y) is an invertible 4x4 matrix depending analytically on is, T (y) and W (v) 
are respectively n x n and (4 — n) x (4 — n) matrices depending analytically on v. 
Additionally 

7» =T + T 1 u + --- , (250) 
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where T has the following Jordan form 

To = Co (In + D ) , (251) 

with /„ being the nxn identity matrix, and Do being a nilpotent matrix, |49| . Section 
6, such that 

D% = 0. (252) 
We would like to show D ^ and, even more, that, 

D 7 *- 1 £ 0. (253) 
Notice that the characteristic equation for T (y) is 

det (T» - C/„) = 0, C = e lfe , (254) 
which, in view of l|251[l . takes the following form 

n-i 

det (T (u) - Ch) = (C - Co)" + E a "~* M (C - Co)" _S + ao M , (255) 

3=1 

where the functions a s (is) for small v have the following expansions a s (v) = a s v + 
(f 2 ) for some number a s , < s < n. Hence, for small v the characteristic equation 
2541) can be recast as 



(C - Co)" + E K-^ + if)] (C - Co)"~ S + ao^ + O (^ 2 ) = 0, (256) 

s=l 

where, importantly, we assume that 

a ^ 0. (257) 

It turns out, that the assumption ao 7^ is equivalent to the following assumption on 
the dispersion relation 

u/") (fco) is finite and nonzero, i.e. < to^ (k ) < 00, (258) 

and the following representation holds 

n\ (iCo)" , x 

To establish this representation we recall that C = e lfc , where k is the wave number, 
and notice that equation l|254(l or l|256|l relate to every v = 10 — loq certain £ ■ (w) and, 
consequently, wave vectors kj (u>) determining dispersion relations. We can also add 
that the algebraic equation 12561) for C = is just another form of the dispersion 
relation (|240|l for u> (k) . Using this observation we can derive H259fl from l|256|) by 
inserting in it C = e lk and Co — e ' fc ° an d, assuming k — kg to be small, we get 

Co [i (k ~ k )] n + a is + O {v 2 ) +0((k- k ) v) = 0, (260) 
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which implies 

a ( w (fc) - u ) = (iC ) n (fc - feo)" + O [(k - k ) n+1 ] . (261) 
Differentiating (|261|1 with respect to k at fc = fco we get 

oowW (fc ) = nl (iC )" , (262) 
implying l|259|) . Notice also that the substitution C = Co in i|255|) yields 

det(T(i/) -Co-fa) = a Q v + O {v 2 ) . (263) 

Recall now that by the Cayley-Hamilton theorem, 015], Section 6.2, any matrix 
T is annulled by its characteristic polynomial, i.e. At (T) = 0. Hence, (|256() holds if 
we substitute C, = T (y) treating all other complex numbers as scalar matrices, i.e. 

n-l 

(T (u) - Co^nf +£ K-^ + (» 2 )] ( T 0) - ( O ln) n ~ s +aovln + {v 2 ) = 0. (264) 

s=l 

Now substituting T [y) = T + T\v + O (v 2 ) into (|264(l and taking in account (|251|l we 
single out the terms linear with respect to v getting the following matrix equation 



(265) 



Suppose now for the sake of argument that (|253|) does not hold, and, hence, Dq ~ l = 0. 
Then in the case of n = 2 we would have Dq = Dq = and the right-hand side of 
the equation (|265|l becomes implying oq — that contradicts the assumption (|257|) . 
Hence, for n = 2, i|253|) holds. In the case of n > 3 the equation (|265|l turns into 



>n— 1 \ ^ r~>n — Srji T^iS— 1 i \ ^ nji-s t- 



(266) 



s=2 



s=2 



so, taking the determinant of the both sides of l|266|) implies, 



det D Co" 1 £ D^T.D^ 2 + £ a^Cr*^" 8 " 1 = (-a ) n 



(267) 



s=2 



s=2 



But, in view of i|252|) . evidently det Do = implying together with (|267() ao — that 
contradicts l|257|) . Therefore, (|253|l is correct and the matrix To = Co C^n + Do) has 
nontrivial Jordan structure. In fact, in view of (|252|l Tq — Co (In + Dq) is similar to 
the Jordan block of rank n, i.e. 



— Co'S'o 



1 1 

'-. 



■•• 











1 



(268) 
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for an invertible n x n matrix Sq. In other words, there exists a basis fo, fi, • ■ ■ , f n 
such that 



Co 1 ^o/o — fo, Co ^b/l — /l +/0i •" j Co ^b/n — /n + /n-l- 



(269) 



The basis fo, fx, ■ ■ ■ , fn reducing To to its canonical form is not unique. What is unique 
is the following set of spans 



Span{/ }, Span{/ ,/i}, Span{/ , /i, . . . , / n -i} • 



(270) 



Possible bases preserving the canonical matrix to the right of So in i|2t)8|l and 
<|2(i9|) are described by the following transformations 



where 



(271) 





" 1 


1 







" 






" 1 


1 







" 







1 















1 








s 















s 


-1 _ 























1 


1 












1 


1 




. 










1 _ 






. 










1 _ 




" 7i 


72 




lu- 


1 


In 

















7i 








ln-1 












s 






















3 


= 1,< 


■ • n 










ll 




72 



























7i . 













(272) 



Let us introduce now a 4 x 4 matrix Q reducing the matrix Q 1 T (0)G (0) to its 
canonical Jordan form Q~ 1 G~ 1 T (0) Q (0) Q. In other words, if we denote 



then we have 



Go V) = G{y)Q 
' T(0) 



Go L (0)T(0)Go(0) = 









W(0) 



(273) 
(274) 



where the both matrices Co 1 T (0) and W (0) are of the canonical Jordan form. Namely, 
T (0) = T takes the following form as in 



Co X T(0) = 



1 1 

'•. 










1 



, and W (0) has canonical Jordan form. (275) 



In the most interesting case of the inflection point for n — 3 the matrix W (0) reduces 
to a scalar. In the case n = 4 there is no W (0), and in the case n = 2 in a generic 
situation W (0) will be just a diagonal matrix. 
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Consequently, the basis fj, j = 0, 1,2, 3 reducing T (0) to the above mentioned 
Jordan form (|274|) . I|275[) can represented as follows 

h = Go (0) bj, j = 0, 1, 2, 3 where (276) 



" 1 " 




" " 




" " 




" " 







1 















, bi = 







1 


, &3 = 




















1 



8. Spectral perturbation theory of the transfer matrix at a point of 
degeneracy 

In this section we develop the spectral perturbation theory for the transfer matrix 
T (lj) defined by (|101|l . This problem has been considered in jHO] for a stationary 
inflection point. For an inflection the essential part of the perturbation theory is 
related to perturbational spectral analysis of the Jordan block of the rank 3, i.e. 



Do = D< 3 > = 



1 
1 




(277) 



Below we extend the spectral constructions from |30) to the case of degenerate points 
of the ranks 4 and 2. It turns out that as in the case of an inflection point, which is a 
degenerate point of the rank 3, the essential part of perturbational spectral analysis 
is reduced analysis of Jordan blocks of the ranks 4 and 2, i.e. 



Do = D 



(4) 



10 

10 

1 





Dn = D, 



(2) 



1 




(278) 



We use the notations 



T (oj) = T (L; oj) ; T (is) = T (w + v) , v=ui- lo . 



(279) 



The transfer matrix T (is) depends analytically on v in a vicinity of v = and it can 
be be reduced to its canonical Jordan form (see pi9" | - j|2"K3 |l . (|2l?8^ and P?l )l -|2"76 ]l ) 



T(is) = g (is) 



T(v) 
W(v) 



Go 1 W 



(280) 



with help of the 4x4 invertible matrix Qq (is) depending analytically on is in a vicinity 
of is = 0. The matrix T (is) in (|280|l is also analytic at is = and has the following 
representation 

T(is) = T Q +T x is + --- , T = Co (I n + D ), (281) 
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where in our case D a — is a Jordan block (|277[1 . (|278[l of the order n — 2,3,4 
correspondingly to the rank of the degenerate point, i.e. 



T(0)=T o = Co (in + D™) - Co 



It is convenient to recast l|281|) as 



1 1 

'■• 



••• 

(n) 











•• 

'•• 1 

1 



T (is) = Co [I* + 1 (i/)] , 5» = + Ti^ + T 2 
T s = Co 2» > s = l,2,.... 

Let us introduce also matrices i^g™' > by 



(2) 





1 



K, 



(3) 







1 



A", 



(4) 









10 



and for every n x n matrix T define the matrix 

T^ = T — [T] nl 2^ n) , where [T] nl is the named entry of T. 
It turns out that the following very special case of T (is) 

1 



To („) = T^ („) = D W + = 



being an exact solution to the equation 

Tp 1 (i/) = vl n 











o ■• 

v 











1 





(282) 



(283) 



(284) 



(285) 



(286) 



(287) 



plays the key role in the spectral analysis of T (is) . For that reason we study first 
spectral properties of To (v). 

Notice that the characteristic equation det (To (is) — (I n ) = for the eigenvalues 
of To (is) is 

C n - is = 0, (288) 

and that the matrix To (is) is a companion matrix of the polynomial C™ — v, |49j. 
Sections 2.2, 2.3. Hence, if we introduce n-th roots of 1 



<^o = 1, ?i 



C2 = s\, 



(289) 
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(290) 



For the most interesting case of an inflection point n — 3 we use another natural 
notation for the roots 

ft = 1, ?i = ?+ = <^ = + ^iV3, ft = *? = ?- = -\ - |i\/3. (291) 

The corresponding eigenvectors of the companion matrix To {v) can be also found, 
and, if one puts them as columns iuanxn matrix &q (v), it takes the form, |5(J|. 
Section I.10-I.13, Section 2.11(Problem 21), 



1 1 


1 


1 1 








2 o 2 


2 2. 







(292) 



Hence, 



T (i/) = ^»6 (i>) A eo 1 (i/) , A = 



1 

?! 








?2 



(293) 



Observe that ©o (^) is a Vandermonde matrix, ^Hjj of order n corresponding to 
n numbers l,<ri,<T2> " • " > an d 



det So (v) = dct 
= (-1) 



11 1 ••• 

1 ?1 ?2 ' ' ' 

1 3 ••• 



n 

1<j <s<n 



(294) 



:» + 2)(^-l) 



-1 ■ (n+2)(n-l)» 



Notice also that 



6, 



1 M = -ej for 1 ) = - 



n 2 v 2 = e 4 n 2 u 2 



1 V n V n 

i -1 -A -2 -2. 

1 q V * <; 1 V n 

i -1 -A -2 -2 



(295) 



?1 ^ r 



1/ n 

-2 S^ 3 - 

?! V » 

-2 2=2 
?2 ^ " 



where a & is conjugate transpose to a matrix S, and £ is the conjugate to a complex 
number £. Evidently the matrix 6 (^) is analytic in z/ 1 /" and 



6 (z/)6 (i/) = i/- V/„. 



(296) 
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Let us consider now the [n, 1] entry of the perturbed matrix T (v) as a new variable 
namely 

5 = I* Mini = E ^ S > ts = P^3i . * > I- (297) 



We will consider the generic case when 



(298) 



The above assumption Ij298(l . as we will show, is equivalent to the fundamental 
assumption i|258|) on the dispersion relation at the point ko- Under the condition 
(|298() the relation (|297|l can be inverted as 



(299) 



where the coefficients t s can be expressed recurrently in terms of tq, q < s. In 
particular 

1 ta 2fi-tifa 



ti = — , t 2 
ti 



t3 



Hence, from (|299|) and (|3UU|) we have 



1 .-. *2 -2 , 2t| — tlt3 _ 3 



v = — v 

ti tf 



(300) 



(301) 



Using the new variable v and (|286|) we recast the perturbed matrix H283JI as a series 



m v: 



T (v) = T 0>)+^z> s T s 



0, s> 1, 



(302) 



To („) = T " } (»/) = DM + ^ = 



1 

o o '■• 



o •• 

v 



•• 

1 





where the matrix %q{v) satisfies also (|292ll . (|293|) . and the matrices T s can be 
expressed recurrently in terms of T^, q < s. In particular, 



T" ~ 
Tr = fi, T 2 
ti 



+ M <f 



(2i| - tit 3 ) tJ 2T 2 T? 



r 1 



In particular, in view of the equalities yield 



n l i 



(303) 



Ti = — , T 2 
Cob 



Cot? + Cot 



2' ^3 



(2t| - tit 3 ) Tf 2T 2 B T, 



jf + fk- (304) 



Based on (|292|l and (|295|l we get the following representation for an arbitrary n x n 
matrix 2t 



9 =0 



(305) 
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where, evidently, (21) are n x n matrices can be found based on the matrix 21 from 
the very relation i|305|) . In particular, one can find that 



Wo = - [sq nl 



t 



t 



Cn-1 Sn-1 



i 



(306) 



f f 


1 


c 2 c 2 




*1 VI 




c 2 c 2 


. c 2 







[2i]„,: 



1 



Cl 



Tn-1 



<Tn-l Sn-lSl 



Sn-lS> 
(307) 

showing as f — + the most significant zero term in the representation (|305J) depends 
only on the entry [2l] nl of the entire matrix 21. This elucidates the special role played 
by the matrix entry [2t] nl . Then from (|293|l . Ij302|l - (|3()7|l we get the following important 
representation 



= A + J2 i)S6 a 1 (P) T S S (v) = A + ^ i 



(308) 



(n-l) 2 
9 =0 



V " 



(n-lf 
s>l g=l 



1/ » 



A + />-(Ti) i +j>" ^Ti /2 



where the matrix Ao is a diagonal matrix defined by (|293[) and, evidently its entries 
1, Si,... , Cn-i defined by l|289[) are all distinct. Notice that the representation l|308[l 
reduces the perturbation analysis of the initial series to the last series in l)308f) . 
The perturbation theory of that series involving a diagonal matrix Ao with different 
elements is much simpler and elementary. The relevant perturbational statements 
needed for the analysis are collected in the following section. 

To analyze perturbations of the matrix Ao we introduce first the following 
auxiliary variable 

Vn = \v, (309) 

Then based on the described results and general facts on the perturbation theory for 
diagonal matrices |46j (the sketch of the theory is presented in the Appendix 2) we 
get 

% (u) = iz>6 (&) e~ s M (Ao + (iz>) Ax + {wf A 2 + • • • ) e^©^ 1 [v) , (310) 
5 (iz>) = (u>) Si + (u>) 2 S 2 + ■ ■ ■ , 
where A s , s > 1 are diagonal matrices. The above formula can be also written in the 
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-S(ii>) 



(Aq + Ax (w) + A 2 (ivf + • • ■ ) e s ^& (i/) , (311) 



form 

- R _1 <3o (i/)e 
5 (iz>) = (i*) 5i + (i/>) 2 5 2 + • • • . 

We would like to point out that in the representation H31UJ) . the eigenvectors 

collected in the Vandermonde matrix &o (v) defined by \292)) are invariant under any 
change of variables described by \2 r Jl\j , \2T0jj . The dependence of the eigenvectors on 
the parameters jj, j = 1, . . . , n comes through terms of proper higher powers of v. 

The representation (|310|l and l|3U9|) imply that the eigenvectors of the matrix T (y) 
(and, hence, in view of \288\) . of the matrix T(v)) are the columns of the following 
matrix 



6 (D)e- s W 



1 1 1 

(iO) ?i (iv) ? 2 (i^) 
(u>) 2 ,l(iv) 2 ^(iz/) 2 



{l n -(iv)S x +0(v 2 )}. (312) 



Observe also that 

detT(i/) = (iz/f detA (l + O(i/)) = (ivf (1 + O (£>)) = i> (l + O (^ V ")) • (313) 
/From [t5S3jl and (f^T5f) we get 

det (T (z/) - C /„) - Co dot 1 (i/) - Co* (l + O (i> 1/n )) . (314) 
Comparing (|314|l with (|263[1 and taking into account (|259|) we get 

n\ (iC ) n 



Co J/ = ao*/, ao 



(315) 



uW (fco)' 

The relation (|315|l combined with H301J) yields the following representation for the 
important quantity ti = [Ti] nl 

n!i" 



ti 



(fco) 



» (fco) 



(316) 



TTie representation yd lb]) , in turn, implies the equivalence of the assumption 
(t\ 7^ 0) to the fundamental assumption \25tfy on the dispersion relation at the point 
fc - Combining (|3H)jl with (|3UT|) and (|3U9j) we get 



W = OiQlV 











, a = 


(fc ) J 



(317) 



The diagonal matrices A s , s > 1 as well the terms of the Taylor series for S (iv) 
can be found recursively (see the Appendix 2 for the details). 

Notice that the eigenvectors of the transfer matrix T (v) in view of H31Q(I and 
JSBJ-CESi) take the form 

©o (v) e- s <^ 
J 4 -„ 



tj (v) = Go (v) 



bi, j =0,1,2,3, 



(318) 



where, according to l|309|l . 

Go iy) = <?o (0) + O (v) = Go (0) + O (v n ) , where n is the degeneracy order. (319) 
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8.1. Spectrum of the transfer matrix at an inflection point 

In this section we derive the asymptotic formulae for the eigenvalues and eigenvectors 
of the transfer matrix T (y) as v = lu — luq — > in the case when the frequency ujq is 
an inflection point, i.e. a degeneracy point of the order 3. We remind the reader that 
in this case according to (J^HU, J2HH), and l(27^l - (|276) we have 



T (v) = Go {v) 



7» 
W{v) 



Go 1 (") . 



(320) 



where Go {y) is a 4 x 4 invertible matrix Go (y) depending analytically on v in a vicinity 
of v = 0, T (y) is a 3 x 3 matrix depending analytically on v in a vicinity of v = 0, 
(f) is a complex valued function analytic in v in a vicinity of v = 0. In addition to 
that, (see H425(l 'l. we have 



T(0) = Co 



1 1 
1 1 
1 



, \W(0)\ = l. 



(321) 



In other words, the basis fj, j = 0, 1, 2, 3 defined by (|276|l reduces T (0) to its canonical 
form 

" Co Co o 
Co Co o 
Co o 
W(0) 



T(0) = S (0) 



a - x (o). 



(322) 



<S. i. J . Eigenvalues of the transfer matrix Observe that it follows from l|293(l . H310|l 
the eigenvalues i] (u), n + (v), r\_ (v) of the matrix I3 + T (v) from H283|) are 



{v) = l + + O (> 3 ) ,r, + (v)=l + c+z> 1/3 + O (> 3 ) , 
r)_ (u) = l + ,-v 1/3 + O (V^ 3 ) , v = t^V/ 3 + O (> 3 ) , 



(323) 



or, in view of l|297|l . 

v- M = i + ti / W /3 + o(i^/ 8 ) ) 

where according to l|291|l 



= 1 



+ t\'\ + v^ + O (> 3 ) , (324) 



1 1. 



Co = 1, ci = c+ = e 1 ^ = -i + -iV3, <r 2 = c 2 = = -~ - ^iV3. 



Notice that we can recast 
Vo 0) = exp 



(325) 



{^ /8 I/V3 + O (>3) j _ ?y+ = exp {tJ^.+ ^/S + O (>3j j ^ (326) 

77_ {v) = exp {t^ 3 ^ 1 / 3 + O (V 2/3 ) } . 
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As follows from the statement (|571f) at least one of r]j (v) must satisfy |ry^ (v)\ = 1. 
Without lost of generality we can choose that one to be rj (y), and, hence for 
sufficiently small 5 > we have 

1% 0)1 = 1 for M (327) 

The representation 13241) together with (|327|l (see also (|3f6|l and (|317[l ) yields 

i 

6 



,1/3 



aoi with a real ao 



V" Oo) 



(328) 



6 



w'" (k ) 



(329) 
(330) 



ti = [Ti] 31 = — api with a real ao 

Hence, l|326|) takes the form 

% {y) = exp {ia ^ 1/3 + O (V /3 ) } , r/+ {y) = exp {ia ^ 1/3 + O (z^ 2/3 ) } , (331) 
r}_ (y) = exp {ia <^ 1/3 + O (V /3 ) } . 

Observe that (|325|l and (|331|l imply that 

if a v > then \r) (y)\ = I, \r) + (v)\<l, \r)_(v)\>l, (332) 
if otov < then \r] (v)\ = 1, \r)_ (v)\ < 1, |r/ + > 1. 



Now, as follows from 
take the form 

If now we denote 



the eigenvalues 6j (v) of the 3x3 transfer matrix T (v) 
j {v)=Q o r lj {u), j = 0,±. (333) 

Co = e iko where fc is real, (334) 



then (EHU-lESU) imply 

O (u) = exp {ifc + k^ 1/3 + O (V /3 ) } = e ifeo+ia ^ I/3 (l + O (V^ 3 )) , (335) 
9+ (u) = exp {ifco + i«o^ 1/3 + O (^ 2/3 ) } - e ik °+ ia °^ /3 (l + O (^ 3 ) ) , 



6- (y) = exp {ifc + ia C-Z/ 1/3 + O (^ 2/3 ) } 
Observe that 



_ ifco+iaoc 



■" 1/a (i + o(, 2 / 3 )). 



if a ^ > then \0 a {v)\ = 1, |<9+ < 1, |0_ {v)\ > 1, 

if a Q v < then \9 Q = 1, |0_ (u)\ < 1, |0 + > 1, 

which can be recast as 

l*o 0)1 = 1, 0)| < 1, \B 0)| > !• 



(336) 



(337) 
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Notice that if Afc = k - fc and Auj = uj - oj = v (|52"g|) . {(TZty and yield 



iAfc = t^V 73 + O (V 2 / 3 ) or Auj = v = -^Afc 3 + O (Afc 4 ) 



implying 



J" (fc ) 

Notice also that (|259|l implies 



: — or iti = a?, — „, ,, , . 
iti w'"(*o) 



6 



It is convenient to introduce 



oj'" (fco) 

1/3 



Co- 



and to rewrite (|335J) and (jSHZjl a s 

O (i/) = exp {ifc +w + (z> 2/3 ) } = e^+™ (l + O (*> 2/3 )) , 
8+ (y) = exp {ifc + + O (z> 2/3 ) } = e ifc °+ k +* (l + O (z> 2/3 
0_ (i/) = exp {ifc +k-v + (z> 2/3 ) } = (l + O (z> 2/3 

The relations 1|342|) . in turn, together with i|32ti[l imply 



(338) 
(339) 

(340) 
(341) 

(342) 



\0o Ml = i, 



- Slg 



n(v) iy)\ > 1. 



(343) 



8.1.2. Eigenvectors of the transfer matrix We recall that we work with the basis in 
which T (0) has its canonical Jordan form as in $281\) . namely 



t(o) = ( 



1 1 
1 1 
1 



Recall also that the eigenvectors e,- (y), j = 0, 1, 2 of the matrix T {y) are the respective 
columns of the matrix So {v) e~ s ( Iiy ) represented by the asymptotic equality (|312() . To 
use (I312|l we need to find the matrix S±. Following the Appendix 2, we first introduce 
a decomposition of a square matrix W into its the diagonal component diag (W) and 
the remaining part W — W — diag (W) with zero diagonal elements, i.e. 

W = [W mj ] = diag (W) + W, where diag (W) = [W m] 5 mj ] (344) 
W = W- diag (W) , 

where 5 m j is the Kronecker symbol. Then we get the following expressions for the 
matrices Ai and Si as follows: 



A x = diag [Si] = 



1 



W m - Wn 



Wi 



n * m; [Si]„„ = 0, (345) 



where W\ = an d Wi = To = 1, w 2 = Si = ?+, w s — <; 2 = 
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Carrying out the operations described in Q345JI . and using (|283|) . I|285[l . (|3l)3|) . I|3U7[) 

we obtain 



1 

3t7 



Pl] 2 l + Pi] 32 
Pi] 21 + ^+ Pl]32 



^+ Pi] 



21 
21 



Pi] 



32 



Pilai + ^+ Pi]32 

?- (pl] 2 l + pl] 32 ) 
Pl] 2 l + pl] 32 



pl] 2 l+<r- pl] 32 

" pl] 32 

Pl] 32 ) 



pl] 2 l 
^+ (pl] 2 l • 



3C ti 



[Ti} 21 + [Ti] 

[Ti] 21 + + Pl] 3 2 



32 



32 



32 

([Ti] 21 + [X 1 ] 32 ) 
^ + [Ti] 21 + [T!] 32 



Fl] 21 

Pl] 2 l 
^+ (Pl] 21 



(346) 

- [rilaa 
H I r i] 32 
Pi] 32 ) 



Ai 



Si 



Pi] 21 + [Tx] E 



3C ti 



1 

3tT 







1 

cr_ 




(347) 



g2pl] 21 +Cl[Tl] 

cipiki+^pi] 



fi-i 




? 2 -i 



l-?2 



T 2 -?l 





f+-l C--1 
T-Ti+<^ + t 2 n ?_-ri+r 2 









-T 2 ?+Ti+T 2 







Pl] 2 l Pl] 32 
' T l = nA , > T 2 = - 



3Coti 



3Coti 



The eigenvectors eo (f), ei (i/) and e 2 (i/) of the matrix %(v) (and, hence, the 
matrix T(is)), corresponding respectively to the eigenvalues ?o = 1> an( l Ci = ?+ = 
i (—1 + i-v/3) and ? 2 = ?- = -H^ 1 ^ 3 ") in vicw of and [|3T7|> take the 

following form 

l + w 2 v + 0(v 2 ) 
e (is) = \v + itiv 2 + O (is ) ,ei(u) = 



1 - ^T 2 V + O (V 2 ) 
±i^+i(T 2 -T 1 )^ + 0(z> 4 ) 

(348) 



y/3-i 



^T 2 i>+0{V 2 ) 

iy^-l , ,'.2 
2 

'.3 



.+ b^=l Tl< )9 + O(0 3 ) 



i^S^ + i( T2 - Tl )^ + 0(^ 4 ) 



e 2 (i/) = 

The equality (|348|l . in turn, implies 
ei (is) - e (is) 



Tl 



P±k 
3Coti 



T 2 



[Ti] 



32 



3C ti 



(*-*) 



T 2 +0 (V) 
l + ^T^ + O^ 2 ) 

^v + 0(v 3 ) 



Tl 



[Ti] 



21 



3Coti 



T2 



U 32 



3Coti 



(349) 



Consequently 



lim eo (is) 
v— >o 



lim 



ei (f ) - e (is) 



6 



T 2 
1 





, T 2 



U 3 2 



3C ti 



(350) 
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Thus we have the following set of eigenvalues and corresponding eigenvectors 



(y) 



(l + 0(z> 2 )), eo(y) = 



6+ {y) = e iko+i ^+ (1 + (z/ 2 )) , e+ {v) = 
0_ (y) = e iko+i ^- (1 + (z> 2 )) , e_ {v) = 



1 + it 2 v + O (v 2 ) 
\v + \t x v 2 + O (ir) 
z/ 2 + i(r 2 - Tl )i> 3 +0(i> 4 ) 

1-'^t 2 v + 0{v 2 ) 

_i+^t_i^±l Tl ^ + ( { ,3) 

'v 2 +i(r 2 -T 1 )^ 3 + 0(z/ 4 ) 

,2\ 
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(351) 



2 



V^5 



1 - i+^r 2 z> + O (z> 2 

i\/3-l. 



2 



3— ■ — r^ 2 + 0(^ 3 ) 



^ 2 +i(r 2 -T 1 )^ 3 + 0(z/ 4 ) 



ri 



U21 



3Co*i 



fT 1 

T 2 = ir^j 1 ; e+ (u) = ei (1/) , e_ (1/) = e 2 (1/) . 



Notice that in view of l|343[l we always have 

\Oo W)\ = 1, l^ig 



(f)| < !> 



(352) 



implying that the vector eo {v) always corresponds to the frozen mode and the 
vector e s jg n (£) (v) always corresponds to the evanescent mode, i.e. the one decaying 
exponentially away from the surface of the photonic crystal. In particular, the two- 
dimensional space Spanjeo {v) , e sign (£) (1/)} describes all possible values of the EM 
field of the ST (scattering theory) eigcnmodes on the surface of the photonic crystal. 
It readily follows from l|35U|) that 



(353) 



Hence, in particular, the two-dimensional space Span {eo (y) , e s j gn (,>) (z/)} , which 
describes all possible values of the EM field of ST eigenmodes on the surface of the 
photonic crystal, converges as v — > to the space Span{/o,/i}, which describes the 
two-dimensional space of all possible values on EM field of ST eigenmodes on the 
surface of the photonic crystal for v = 0, i.e. at the frequency luq of the frozen 
eigenmode. 

Hence, in view of l|318|l and (1319(1 . we have the following representation for the 
eigenvectors tj (v) = Go (u) ej (v), j = 0, 1, 2, 4 of the transfer matrix T (v) 





" 1 " 




' " 


Span{/ ,/i}, /o = 





, h = 


1 













eo (f) = Go (0) 



1 + \T 2 V + O (V 2 ) 

\v + ir^ 2 + O (z/ 3 ) 
-z> 2 + O (O 3 ) 



0(v 6 



U21 



d (1/) - e+ (1/) - Go (0) 



i+V3 
2 



1 - 

i+V3 ,<, _ h/3+1- .-.2 



3Coti 

T 2 Z> + O (Z/ 2 ) 



^2 



U32 



3C ti 



(354) 



2 -v — 1 1 " ti^ 2 + O (z> 3 ) 
^iiz> 2 + O (z> 3 ) 
" 0(z> 3 ) 
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e 2 0) = So (0) 



v/3 



T 2 V + O (v 2 



l=g& + o (u 3 ) 
' 0(y*) 



c 3 (y) = g (0) 



1 + O (*> 3 ) 



(355) 

Combining now (|354ll . (|355l) with (|319[l and (|276[) we get the following representations 
for the eigenvectors tj (v) 

to 0) = (1 + it 2 v + O (v 2 )) f Q + (w + \t x v 2 ) fi - v 2 \ 2 + O (v 3 ) , (356) 
e+ (i/) = ei (i/) = (l-^^^^ + O^ 2 ) ) f + 



i + \/3, iV3 + l 



r^ 2 fi + 



iV3 + l 



^ 2 f 2 + O (*> 3 ) 



M = e 2 (f) = ( 1 + + O (u 2 )) fo + 



\/3-i, i-v/3-1 



TlV 2 fl + 



or, using the symbols <;± = — | ± from (|291[1 we have 
e (v) - (1 + ir 2 j/ + O (z> 2 )) fo + (h> + iri/> 2 ) fi - *> 2 f 2 + O (v 3 ) , 



(357) 



e+ (v) = d (y) = (1 + k+r 2 j/ + O (;> 2 )) f + (k + u + st^ 2 ) - <^_z/ 2 f 2 + O (;> 3 ) , 
e_ (i/) = e 2 (i/) = (1 + k_T 2 z> + O (v 2 )) f + (fc-0 + c+tiv 2 ) fi - <; + ^ 2 f 2 + O (f 3 ) . 



Notice that the equalities (|357|l imply 

lim e (u) = lim e± {v) = fo, 



(358) 



indicating, in particular, that the three vectors eo (y) and e± (v) become almost parallel 
as v — > 0. To have a nicer way to trace the two-dimensional spaces Span {eo iy) , e± iy)\ 
we introduce the following two vectors 



(?+ - 1) 

e-(,)-e (.) = 
if (?_ - 1) 

Notice that the equalities (|354|l . (|359|l imply 



C. — i 



fi + :/ X ^^ f 2 + O (^ 2 ) ; 



fi + 



lim fct H=r 2 f + f 1 , r 2 = Pf . 
i/— o 3C ti 

Then the relations (JSSHJ), (|335)l and yield 

Span{e (f ) , e± [v)} = Span{e iy) , f)± (v)}. 
lim Span {e (f ) , e± [y)} = Span {f , fi} . 

v— >0 



- 1) 



(359) 
*fa + O {v 2 ) 

(360) 
(361) 
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8. 2. Spectrum of the transfer matrix at a degeneracy point of the order 4 

In this section we derive the asymptotic formulae for the eigenvalues and eigenvectors 
of the transfer matrix T (v) as v = u) — wq — ► in the case when the frequency coq 
is a degenerate point of order 4. In this case since n — 4 the transfer matrix T {v) 
defined by l|242JI is such that T (0) is a Jordan block of order 4, as it follows from the 
analysis carried out in the previous Section. The mentioned analysis implies also that 
the matrix T (v) defined by (j249j) is a 4 x 4 matrix and that T (v) = T (v). 

Without loss of generality we assume that loq = lj (kg) is a point of local minimum 
of the dispersion relation u> (k) in a vicinity of kg. In the later case for v — uj — ujq > 
and v small there must be two propagating Bloch modes and two evanescent modes. 
Consequently, there will be two eigenvalues of the matrix T (y) with absolute value 
1, one eigenvalue with absolute values lesser than 1 and one eigenvalue with absolute 
value larger than 1. 

8.2.1. Eigenvalues of the transfer matrix Observe that it follows from Ij293(l . H^lOjl 
that the eigenvalues T]a (l>), j — 0, 1, 2, 3 of the matrix I4 + 1 (v) from l|283|) are 

r,j (u) = 1 + + O (V /2 ) , J = 0, 1, 2, 3, (362) 

or, in view of l|297jl . 

Vj M = 1 + ^%v 1/4 + O (v 1/2 ) , i = 0, 1, 2, 3, (363) 

where 

?j = j = 0, 1, 2, 3. (364) 

Notice, we can recast (|363ll as 

Vj (y) = exp {tj 7 V 1/4 + O (^> 1/2 ) } , j = 0, 1, 2, 3. (365) 

As we have found at the beginning of the section, \r]j = 1 for exactly two values 
of j = 0, 1, 2, 3. Using and (|3T7jl we get 

> 0, (366) 
and 

> 0. (367) 
z^ 4 = iaov 1 / 4 + O (v 1 ' 2 ) . (368) 
Hence, l|365|) takes the form 

{v) = exp {ia ^^ 1/4 + O (i> 1/2 ) } , j = 0, 1, 2, 3. (369) 
Observe that for ao > 0, l|369() implies that 



aoi with a real ao 



24 

T^Tko) 



[1 



U41 



a with a real ao 



24 



u> 



W (fco) 



l%WI = l%WI = i, l%H<i> 1% (")!>!• 



(370) 
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Now as it follows from 
take the form 

If now we denote 

then $£$-$22$ imply 



the eigenvalues 6j (v) of the 4x4 transfer matrix T (y) 



9,>) = Co^>), 3 =0,1, 2, 3. 
Co = elk ° where fc is real, 



6j (u) = exp {ifco + k*o^ 1/4 + O (V /2 ) } 

jk„ +iao ^ f 1 + Q ( v i/2\ \ j = 0, 1, 2, 3. 



Observe that 

\6 Ml = |*2 Ml = 1, |fliMI<l, |6>3 (^)l > I- 
Notice also that (|259H implies 

24 



4a4 

a = "oCo = 



;( 4 > (Ad) 



It is convenient to introduce 



and to rewrite (|373|l as 



1/4 



(371) 
(372) 

(373) 

(374) 

(375) 
(376) 



iy) = exp {ifc + kji> + O (v 2 ) } 



(l + 0(*> 2 )), j = 0,1,2,3. (377) 



8.2.2. Eigenvectors of the transfer matrix We begin this section with the reminder 
that we work with the basis in which T (0) has its canonical Jordan form as in y281)) . 
namely 

"110 
110 
11 
1 



r(o) = c 



Based on (|2"%5|) . i|3TUJ> . l|3T7|) and lj3T%|) we can find the eigenvectors (i/) of 

T (v) corresponding to its eigenvalues 9j (v). 

Recall also that the eigenvectors ej (y), j — 0, 1, 2, 4 of the matrix T (v) are the 
respective columns of the matrix ©o (y) e~ s ^ represented by the asymptotic equality 
<|312|) . To use <|312[) we need to find the matrix Si. Following the Appendix 2, first we 
introduce a decomposition of a square matrix W into its diagonal part diag (W) and 
the remaining part W = W — diag (W) with zero diagonal elements (as in H344|) ) 



W = [W„ 



diag (W) + W, where diag (W) = [W mj S mj ] , W = W - diag (W) . 



where 8 m j is the Kronecker symbol. Then we get the following expressions for the 
matrices Ai and Si 



Ai = diag(Wi), [5i] nr 



1 



Sm— 1 Sra — 1 



n 7^ m; [Si 



0, (378) 



where W\ = (^1/ »?o = 1, ?i = i, ?2 = -1, ?3 = -i- 
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Carrying out the operations described in (|37Sj) . and using <!303|) - P77jl 

we obtain 



Ti + r 2 
?i (ri?i + t 2 ) 

<^2 (Tlf2 + T 2) 
<^3 (Tl?3 + r 2) 



Pi] 



31 



4C ti 



^2 



[Ti] 



Ti + r 2 ?i 
«i(ti+t 3 ) 

<^2 (Tl?2 + T2?l) 
?3 (TlC3 + T 2<?l) 



42 



4C ti 



7" 1 + T 2<?2 Ti + T 2<?3 

t 2 <; 2 ) ci (n?i + r 2 <j 3 ) 

T"2) <?2 (ri? 2 + T 2<^) 
T"2?2) <^§ (Ti + 1"2) 



?i (nci 

<^3 (Tl£3 



(379) 



5i = 





Ci(Ti^i+r 2 ) 

l-f 2 
C3(Tl^3+T 2 ) 



Ti+T2fi 



< ;2(ti*;2+t2';i) 

, fl _ f2 
f3(Tlf3+T2f l) 



52-1 
£l(211il+l2£2) 



£3_(211i3+l2£2) 
<;2-?3 



f2(Tlf2+ , ''2f 3) 
?3-l 





(380) 



1— fa fi-f3 

The eigenvectors ej {y), j = 0, 1, 2, 3 of the matrix T (f) (and, hence, the matrix 



T{y)), corresponding to respectively the eigenvalues 
view of H312|) and l|347|) take the following form 



exp 



(2fl), j =0,1,2,3 in 



eo (f ) = 



e 2 (1/) 



l + i(ri +3t 2 )^ + 0(;/ 2 ) 
i£ + (n - t 2 ) i/ 2 + O (v 3 ) 
-;> 2 + i (3n + t 2 ) ^ 3 + O {^) 
-w 3 + 3 (n - t 2 ) z> 4 + O (;/ 5 ) 

1 - i (r x + 3r 2 ) z> + O (z> 2 ) 
-i^+(r 1 -T 2 )i/ 2 + 0(r> 3 ) 
-;> 2 - i (3n + t 2 ) z> 3 + O (*> 4 ) 
iz> 3 + 3(r 1 -r 2 )i/ 4 + 0(z/ 5 ) 



Tl 



[Ti] 



31 



r 2 



U42 



SCoti 8C ti 
The equality Q381fl . in turn, implies 



ei (1/) 



63 M 



l-(r 1 +3r 2 )z/ + 0( ! / 2 ) 
-/^+(T 2 -Ti)^ 2 + 0(r> 3 ) 
£ 2 + (3ti+t 2 )i> 3 +0(z> 4 ) 
-z> 3 + 3 (n - r 2 ) z> 4 + O (z> 6 ) 
(381) 

l + (r 1 +3r 2 )// + 0(;> 2 ) 
^+(r 2 -r 1 )// 2 + 0(^ 3 ) 
z> 2 - (3ti + r 2 ) z/ 3 + O (i> 4 ) 
£ 3 + 3 (n - r 2 ) i> 4 + O (£ 5 ) 



ep (g) ~ gx (y) 
(l + i)z> 



Consequently, 



n + 3t 2 + o (0) 

l + (l-i) (r 1 -r 2 );/ + 0(;> 2 ) 
(i - 1) £ + (3ti + r 2 ) z/ 2 + O (v 3 ) 
-\{? + O (v 3 ) 



(382) 



lim eo (f) 

t/— >o 



1 





SCoti 



lim ej ( ; } - e ° ( " } 



Tl + 3t 2 

1 






T 2 



8C ti 



3 = 1,2,3, 



(383) 
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In view of the above and i|377|) , we have the following complete set of eigenvalues 
and corresponding eigenvectors of T (u) 



9 (is) = e ik ° +i »(l + O(i> 2 )), eo(i/) 



9 1 (u) 



ri 



o ifc +ii><;2 



(l + 0(z> 2 )), ei (i/) = 



(l + 0(z> 2 )), e 2 H 



(l + 0(z> 2 )), e 3 (i/) 



l + i(rx + 3r 2 )!/ + 0(z/ 2 ) 
i£ + (n - t 2 ) z/ 2 + O (^ 3 ) 

-z> 2 + i (3n + t 2 ) v Z + O (z/ 4 ) 
-iz> 3 + 3 (ri - r a ) ^ 4 + O (z> 5 ) 
l - (n + 3t 2 ) £ + O (^ 2 ) 

-v+ (r 2 -r 1 )^ 2 + 0(^ 3 ) 
z> 2 + (3n + r 2 ) z/ 3 + O \v 4 ) 
_ -i> 3 + 3(r 1 -r 2 )^ 4 + 0(^ 5 ) 

1 - i (n + 3t 2 ) v + O (z> 2 ) 
-\v + (n - t 2 ) z> 2 + O (z> 3 ) 
-z> 2 - i (3n + r a ) *> 3 + O (z> 4 
iz> 3 + 3 (n - r 2 ) i> 4 + O (v 5 ) 

I + ( Tl + 3t 2 ) v + O (i> 2 ) 
v + (r 2 - ri) v 2 + O (z> a ) 
z> 2 - (3n + r 2 ) z> 3 + O (z> 4 ) 
£ 3 + 3 (n - r 2 ) £ 4 + O (z> 5 ) 



(384) 



."■1 



T 2 



1J42 



8C ti 8C ti 
Notice that in view of l|374|l we always have 

\e Q {v)\ = i, \e x {v)\ < 1, 



(385) 



implying that the vector eo (v) always corresponds to the frozen mode and the vector 
e\ {v) always corresponds to the evanescent mode, i.e. the one decaying exponentially 
away from the surface of the photonic crystal. 

Hence, in view of i|318[l and (|319() . we have the following representation for the 
eigenvectors tj (y) = Go {v) &j {y), 3 — 0, 1, 2, 4 of the transfer matrix T [y] 



eo M = Go (0) 



ei (1/) = Go (0) 



1 + i (n + 3r 2 ) v + O (v 2 ) 
w+{t 1 - r 2 ) v 2 + O (v 3 ) 
-v 2 +i(3r 1 +r 2 )^ 3 + 0(p a ) 



-\v 3 + O (v 4 ) 



(386) 



1 - ( Tl + 3t 2 ) v + O (v 2 ) 
-v+ (r 2 - ri)^ 2 + (z/ 3 ) 



! > 2 + (3r 1 +T 2 )^ 3 + 0( i / 4 ) 

'.3 , /-> / >.4\ 



-z> 3 + o(z> 4 ) 
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e 2 0) = Go (0) 

e 3 (i>) = So (0) 



l-i(Ti + 3T 2 )u + 0(i> 2 

—10 + (Tl - T 2 ) Z/ 2 + O 

-^ 2 - i (3n + t 2 ) z> 3 + o {0 4 ) 
i0 3 + O (0 4 ) 

1 + (n + 3r 2 ) + (0 2 ) 
0+{t 2 ~t 1 )0 2 + O{0 3 ) 
2 - (3n + t 2 ) 3 + O (0 4 ) 
3 + O (0 4 ) 
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(387) 



Combining now (3gg) , (p?8Tj) with fiTfy and (|27H|l we get the following 
representations for the eigenvectors tj (v) 

to {v) = [1 + i (n + 3r 2 ) + O (0 2 )] f + [i0 + (ti - r 2 ) z> 2 + O (;> 3 )] f x + (388) 

[-0 2 + i (3ri + T 2 ) 3 ] f 2 - iz> 3 f 3 + O (0 4 ) , 
ei (f) = [1 - (T! + 3r 2 ) + O (^ 2 )] fo + [-0 + (t 2 - r x ) z> 2 + O (*> 3 )] fi + 

[;> 2 + (3n + r 2 ) z> 3 ] f 2 - z> 3 f 3 + O (;> 4 ) , 
e 2 (i/) = [1 - i (n + 3r 2 ) // + O (^ 2 )] f + [-i* + (n - r 2 ) */ 2 + O (0 3 )] h + 

[-0 2 - i (3n + r a ) *> 3 ] fa + i^ 3 fs + O (0 4 ) , 
e 3 (i/) = [1 + (n + 3r 2 ) z> + O (^ 2 )] fo + + (r 2 - r x ) z> 2 + O {0 3 )\ h+ 

[0 2 - (3ti + r 2 ) ^ 3 ] fa + /> 3 f 3 + O (;> 4 ) . 

It readily follows from Ij384|l and Ij388(l that 



lim e,- fi/) 

j/->0 J 



>o it/ (<;j — 1) 



, J -0,1,2,3, 



ti + 3r 2 
1 





.7 = 1,2,3. 



(389) 



lim n e i O) = fo, i = 0,1,2,3, 

v — >0 



iim - e ° ( ; } 

v— >o if (<^j — 1) 



(390) 



(ri+3r 2 )f +fi, j = 1,2,3, 



indicating, in particular, that all four vectors tj (y), j — 0,1,2,3 become almost 
parallel as v — ► 0. To have a nicer way to trace the two-dimensional spaces 
Spanjeo (v) , ei (y)} we introduce the following vector 



e (t/) - ei (y) 
0(i + l) 



ei (zj - e (v) _ ti (v) - e (v) 
10 (<Ti - 1) iz> (i - 1) 

= (ri + 3r 2 ) f + fi + O {0) . 



(391) 
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Notice that the equalities (|390|l imply 

limh 1 ( i /) = (T 1 + 3r 2 )f + f 1 . (392) 

u—>0 

Then the relations (|3^U|l - l|3l?2l yield 

Span {e {v) , e x (v)} = Span {e (v) , f)i (v)} , (393) 
lim Span{e {v) , e x {v)} = lim Span{e (^) , t)i {v)} = Span{f ,fi} . 

8. 3. Spectrum of the transfer matrix at a degenerate point of order 2 

In this section we derive the asymptotic formulae for the eigenvalues and eigenvectors 
of the transfer matrix T (v) as v = u) — ujq — > in the case when the frequency ujq is a 
degenerate point of order 2. In this case since n = 2 the transfer matrix T (v) defined 
by l|242|) is such that T (0) is a Jordan block of 2, that follows from the analysis carried 
out in the previous Section. 

Without loss of generality we assume that u>o = w (ko) is a point, say, of local 
minimum of the dispersion relation to (k) in a vicinity of ko. In the later case for 
v = lu — ujq > and v small there must be two propagating Bloch modes associated 
with the chosen branch to (k) and two modes which either both are propagating or 
both are evanescent modes. 

8.3.1. Eigenvalues of the transfer matrix Observe that it follows from 1|293|) . i|310[l 
the eigenvalues r/ (v), n 1 (v) of the matrix I2 + 1 (v) from Ij283(l are 

r) (p) = 1 + Vt + 0{i>), ri 1 {v) = l-y/Z + 0(u), (394) 

or, in view of l|297() . 

r) (v)=l + VhV^+O(v), r 1l {v) = l-^/i x ^+0{v), (395) 

where, we recall, according to l|291|l 

Co = 1, ft = -1- (396) 

Notice, that we can recast (|395H as 

=exp{Vti"V^+0(^)}, (1/) = exp{-^y/v + 0{y)}. (397) 

Since we have two propagating modes for sufficiently small positive v there exists a 
sufficiently small S > such that 

1% MI=»? 1 M = 1 for |i/|< 5. (398) 
Using lEHjJl and i|3T7|) we get 



y 7 ^" = aoi with a real a — , ^ > 0, (399) 

V w ( fc o) 



ti = 



[Ti] 31 = — a§ with a real ao = 
iceo^/v + O {u) 



(400) 
(401) 
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Hence, (|397|) takes the form 

7y (V) = exp{ia V^+ O (v)} , r) x (y) = exp {-\a \/v + O (u)} . (402) 

Now as follows from (|283(1 . the eigenvalues 9j (v) of the 2x2 transfer matrix T (y) 
take the form 

*j M - Co»?i M . j = o, i. 

Recalling that 



S.0 — e ' 

where fco is real, we get from (|4U2|I - (|4U4|) 

00 (y) = exp {ifc + ia ^ + O (i/)} = e ife °+ iQ! °v^ (1 + (*,)) , 

01 (i/) = exp {ifc - fo ^ + O (i/)} = e ifc »- ia »^ (i + (z/)) . 

Notice also that (|259|l implies 

2 



(403) 
(404) 



(405) 



«o 



-"oCo 



w" (fc 



It is convenient to introduce 

and to rewrite (|335|l and (|337|l as 

9 (y) = exp {ifc +w + (v 2 ) } = e ik « +iC/ (l + O (v 2 )) , 
0i (i/) = exp {ifc - fe> + O (i^ 2 )} = e ifco_u> (l + O (^ 2 )) . 



(406) 
(407) 
(408) 



8.3.2. Eigenvectors of the transfer matrix. We recall that we work with the basis in 
which T (0) has its canonical Jordan form as in y281)) . namely 



T(0) = Co 



1 1 
1 



Based on (128311 . I|310|l . (|311|l . I|317|) and l|318|) we can find the eigenvectors ej {v) of 
T (y) corresponding to its eigenvalues 6j (y). 

Recall also that the eigenvectors e.j{v), j = 0,1 of the matrix T (v) are 
the respective columns of the matrix (So (v) e~ s ^ represented by the asymptotic 
equality l|312|) . To use (|312|) we need to find the matrix Si following the section 
"Perturbation theory for diagonal matrix" in the Appendix. Namely, first we introduce 
a decomposition of a square matrix W into its the diagonal component diag (W) and 
the remaining part W = W — diag (W) with zero diagonal elements, i.e. 



Ai = diag(Wi), [Si 



-1 — 



Wi 



n * m; [Si] nn = 0, (409) 



where Wi = , = 1, ?i = -1, 



Carrying out the operations described in (|409|l . and using (|283|) . I|285|l . (|303l) . 1307( 1 

we obtain 



Ti 



U22 



2C ti 



1 -1 
-1 1 



Si 



Fl]22 

4Coti 



1 

-1 



(410) 
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The first two eigenvectors eo (v) and e\ (v) of the matrix T (y) (and, hence, the matrix 
T (v)), corresponding to eigenvalues <r = 1 and <;i = i in view of (|312fl and i|347fl take 
the following form 



eo {v) 



l + Wv + 0(v 2 ) 
\v + \ti? + O (v 3 ) 

Pi] 22 
4C ti ' 



l-h-v + O (v 2 ) 

-IV + \TV 2 + O (l> 3 ) 



(411) 



Thus using l|411|) and (|405|l we get the following complete set of the eigenvalues and 
corresponding eigenvectors of T (u) 



9 (v) = j k ° +i *(l + O(v 2 )), eo(u) 
9 1 (v) = e ik °- i *(l + 0(v 2 )), e x {v) 



1 + itv + O (v 

\V + YTV 



2 ' 0(v 



(412) 



l-h-v + O {v 2 ) 
-\v + iru 2 + O (v 3 ) 



Hence, in view of i|318|l and (|319|) . we have the following representation for the 
eigenvectors tj (v) = Qq (v) ej (y), j — 0, 1, 2, 3 of the transfer matrix T (v) 



eo M = Go (0) 



1 + \tv + O (y 2 
w + O (v 2 ) 



, ei (v) = Go (0) 



e 2 (v) = Go (0) 



0(v 2 ) 
0(v 2 ) 
1 + O (v 2 ) 



O £>' 



, e 3 (y) = Go (0) 



1-iTV+O {y 2 
-w + 0(v 2 ) 
Q(v 2 ) 
0(v 2 ) 

0(v 2 ) 
1 + O (v 2 ) 



(413) 



(414) 



Combining now 14131) with l|319|) and (|276[1 wc get the following representations for 
the eigenvectors tj {v) 

e {v) = [1 + \tv\ f + i^fi + O {v 2 ) , ei {v) = [1 - irz>] f - i*f a + O (z> 2 ) , (415) 
e 2 (v) = f 2 + O (z> 2 ) , e 3 (v) =h + (v 2 ) . 

It readily follows from l|413|l and l|415|) that, 



lim e 3 (y) 

v^O 



, j = 0, 1; lim 

jy^O 



ei (z/) - e (z/) 
-2i£ 



[Ti] 



22 



4C ti 
Pl] 22 



lim ej (z/) = fo, J = 0, 1; lim — = rf + fx, r = 

!/^0 1/^0 —211/ 4(,gtl 



(416) 



(417) 



indicating, in particular, that the two vectors eo (f) and ei (v) become almost parallel 
as v — ► 0. To have a nicer way to trace the two-dimensional spaces Span {eo (y) , ei (v)} 
we introduce the following vector 



ei (i/) - e (y) 
-2iz> 



(418) 
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Notice that the equalities (|417|l and (|418|l imply 

Iim hi(i/)=rf +fi, r=^f . (419) 

Then the relations (|4THll - ll4T^) yield 

Span {e (u) , d (i/)} = Span {e (f) , f)i (i/)} , (420) 
lim Span{e (v) , e x (i/)} = lim Span{e (i/) , f)i (i/)} = Span {f , fx} . 

9. Transfer matrix and the flux 

9.1. Transfer matrix and the flux for an inflection point 

The case of an inflection point is the case with degeneracy index n = 3. In this case 
according to (|249[1 . the matrix T (y) is a 3 x 3 matrix and W (y) is a 1 x 1 matrix, i.e. 
just a complex number W (y). Let the matrix T (y) be defined by l|249[l . Consider 
first the matrix at the frequency of the frozen mode ujq, i.e. for v = 0: 

T (0) = T = Co {h + D Q ) , |Col = 1, Dl = and D 2 £ 0. (421) 

Then there exists a canonical basis /o, /i, f% related to the matrix D such that 



£>o7 2 = 0, fa = DSh, /i=A>/ 2 . 



(422) 



In fact, the basis /o, fi, fi is not unique and is defined up to some transformations. 
The equalities (|421|) and l|422|) imply the following representation for To 



Co 1 ^o/o — /oj Co ^b/i — /i + /oj Co ^0/2 — /2 + /l- 
or, in the basis /o, A, we have 



(423) 



7o = Co 



1 1 
1 1 
1 



detT = Co- 



(424) 



In addition to that, (|2l3f and |gZH imply detT(» = detT(0) IT (0) = CqV (0) = 1, 
and, since |Col = 1 we get 



IT(0) = Co 3 , |W (0)| = 1. 



(425) 



We reiterate that our fundamental assumption is that Co is triply degenerate and 
IT (0) is an eigenvalue of T (0) that differs from Co, i-e. 



^(0)^Co- 



(426) 



Observe now that since according to g55j) |W(0)| = 1 then IT (0) = [IT (0)] 1 and, 
hence, the relation H426JI can be rewritten as 



Co IT (0)^1. 



(427) 
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Recalling again the relation (|249|) between the 3x3 matrix T (0) and the original 
4x4 matrix T (0) we introduce a basis f , f i , fa , f3 in the four-dimensional space such 
that 



Si = P3S- 1 (0) fj, j = 0, 1, 2, where P 3 



X 2 

x 3 
x 4 



Xi 

x 2 

X; 



(428) 



and the fourth vector f 3 is the eigenvector T (0) related to the eigenvalue VF(0). 
Evidently the vectors fo, fi, fa are 4-dimensional representation of respective vectors 
/o, fi, /a- Based on the above and the relations i|423l) we get 



^ofo — Cofo, %h — Cofi + Cofo, %h — Cofz + Cohi 



(429) 



Tof 3 = W^(0) f 3 . 



(430) 



In particular, the matrix 7o has just two genuine eigenvectors fo and f 3 with 
corresponding distinct eigenvalues Co aim W(0). The vectors fo and f 3 correspond 
respectively to the frozen mode and the only propagating mode at the origin. The 
vectors f 1 and f 2 correspond respectively to linearly and quadratically growing Floquet 
modes at the origin. 

Recall now that Tq is a J-unitary matrix, i.e. 



T$J% = J, [To$i,T $2] = [$i,$ 2 ] for any $x,$ 2 . 
Using 142911 and l|431|l we obtain the following identities 

[fo, fi] = [T fo, Toft] = [fo, fi + fo] = [fo, fi] + [fo, fo] , 



(431) 



(432) 



[fo, fa] = Pofo, T fa] = [fo, fa + fi] = [fo, fa] + [fo, fi] , 
readily implying that 

[fo, fo] =0, [f ,fi] =0. 
Then using g2§J|, P3T|) again we get 

[fx, fa] = [Tofi.Tofa] - [fi + fo, fa + fi] 

= [fi,f 2 ] + [fo,f 2 ] + [fi,fi] + [fo,fi]. 

The equalities igHH together with yield 

[fo,f 2 ] + [fi,fi] =0. 

and, since [fi, fi] is a real number, we consequently have 

[fo,f 2 ] = [f 2 ,fo] =-[fi,fi], lm{[f ,f 2 ]} = 0. 

We also have the relation 

[fi.fi] = [TohM = [fi+fo,fi+fo] 

= [fi,fi] + [fo,fi] + [fi,fo] + [fo,fo], 



(433) 
(434) 

(435) 

(436) 
(437) 
(438) 
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but, in view of (|434p . it is already satisfied and does not produce a new relation. The 
remaining relation is 



[f 2) f2] = [T f 2) Tof 2 ] = [f 2 + f 1 ,f 2 + f 1 ] 

= [f 2 ,fe] + [fi,fa] + [fa.fi] + [fi.fi] 



yielding 



[fi,fa] + [fa,fi] + [fi,fi]=0. 
Notice that for a natural number m > 2 we have 



(439) 



(440) 



" f 


1 


" 


m 


" t 


m 


(m— l)m 
2 







1 


f 







1 


m 


(441) 








f 










1 





This identity together Ij429|) and l|431|) imply 

[Co 1 ^o] m fo = fo, Vh^h+mfo, 

r i _ 1 _ r m, , , (m— l)m, 
[Co To] f 2 = f 2 + mfi + ^ 5^— fo, 



([Co 1 ro] ro ) t J[Co 1 ro] m = J. 



(442) 



(443) 



Notice that the vectors 7^™fi and T™fe representing the EM field at points mL grow 
respectively linearly and quadratically asm^oo. 
Using l|442|) and Ij443|l we obtain 



[fa, fa 



(m — l)m, , , (m-l)ra, 
fa + mfi + * fo, fa + mfi + - s-^— fo 



= [fa, fa] + m 2 [fi, fi] + m [f 1; f 2 ] + m [f 2 , h 



(444) 



(m — 1) m 



([fo,fa] + [fa,fo]), 



implying 



(m — 1) 

m [h,h] + [fi, fa] + [fa, fi] + ^r-^ Qfo, fa] + [fa, fo]) - 0. 



Combining l|445|l with l|437|l we get 

[fi,fi] + [fi,fa] + [fa,fi] =0. 
which is identical to (|440|l . Consider now 



[fa, fi 



implying 



(m — 1) m 
f 2 + mfi H f , fi + mf 

[fa, fi]+m[fi,fi] +m[f 2 ,fo], 
[fi,fi] + [fa,fo] =0, 



(445) 



(446) 



(447) 



(448) 
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which is equivalent to (|436[1 . So, consideration of powers 7^ m of the transfer matrix 
have not produced new identities. Observe now that H427|) . I430|) and l|131|) imply 



[fo,fs]=0. 

Collecting P34) . lEE37?|) . lEHU|) and lEHSI) we get the following system 



[f 1 , f 2 





[fo, fo] 


= 




[fo, fl] 


= 


[fo.faH 


- [fl,fl] 


= 


[f 2 ,fi]4 


- [fl, fl] 


= 




[fo, [3] 


= 



Notice that since J is Hermitian l|452|) implies 

Im{[fo,f2]} = Im{[f 1 ,f 1 ]} = ! [fo,f2]=Re{[f ,f 2 ]}. 
In addition to that, l|452|) and l|453j) yield 

[fo,f 2 ] =Re{[f ,f 2 ]} = -[fi,fi], 



(449) 



(450) 
(451) 
(452) 
(453) 
(454) 



(455) 



(456) 



2Re{[fi,fa]}=-[fi,fi]. (457) 

Let us show now that 

[fi,fl] ^ 0. (458) 

Indeed, assume for the sake of the argument that [fx, fi] = 0. Then, in view of l|450|) - 
(|452[) and l|454l) we have 

[fo.fi] =0, i = 0,l,2,3, (459) 

or, in other words, 

(Jfo.fi), 3 = 0,1,2,3. (460) 

Since fo, fi, f 2 , f3 is a basis in the 4-dimensional space the relations l|460|) imply that 
Jfo = 0, and, consequently, that fo = since evidently J is an invertible matrix, 
fo = is impossible, and we must conclude that the relation Ij458(l holds. 

Since [fi,fi] is the flux corresponding to the Floquet mode described by f\ the 
relation \458j) signifies a fundamental fact that the Floquet mode described by f\ has 
nonzero flux. 

Observe that (|456|l . (|457|l and (|458|l imply 

2 Re { [f x , f 2 ] } = Re { [f , f a ] } = - [f 1 , f 1] ± 0. (461) 
Notice also that in view of 1450|l . 145 1|) we have 

[afo + Pfi,afo + Pfx] = \P\ 2 [fufi] , (462) 



[/a + af o +0fi,f 2 + afo + 0h] (463) 
= [f 2 , f 2 ] - 2 Re {a} [h,h] + 2Re{/3 [/ 2 , /J} + |/3| 2 [/i,/i] , 
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and, hence, we have: 

\u u\ 

as u runs over Span(/o,/i) ' , then runs over [0, +00) , (464) 

L/i)/iJ 

as u runs over Span (/o, fi, $2) r ' , , then runs over (—00, +00) . (465) 

[/1, Ji] 

The relation (|465|l follows from i|463|) if we set (3 — and let a run over all real values 
(—00, +00). In other words, for all vectors u from the Span (fo, f\) the corresponding 
fluxes have the same sign, whereas in the case of Span (/o, /1, $2) the flux can be any 
real number. 

9. 2. Transfer matrix and the fluxes for a degenerate point of order 4 

In the case of a degenerate point of order 4 the transfer matrix becomes a Jordan 
block of rank 4 and according to l|269(l there exists a basis fj, j = 0, 1, 2, 3 in C 4 for 
which we have 

r f = Cofo, %h = Cofi + Cofo, %h = Cof 2 + Cofi. ( 466 ) 



%h = Cofs + Cfa. (467) 

Notice that the three equations Q466J1 are exactly the same as the three equations l|429|l 
for the inflection point. Hence, the identities l|450|) - (|453|) in this case are 

[fo.fo] =0, (468) 

[fo,fi]= 0, (469) 

[fo,f 2 ] + [fi,fi]=0, (470) 

[fi,f 2 ] + [f 2 ,fi] + [fi,fi] =0. (471) 

Using now (|431(l . 1466^ and (|467(l we obtain 

[fo, fs] = [Tofo, T f 3 ] = [fo, fa + h] , (472) 

implying 

[fo,fa]=0, (473) 
which together with H470jl and (|471fl yields 

[fi,fi] = 0, [fi,f 2 ] + [f 2 ,fi] = 0. (474) 

Observe an important difference of the case of a degenerate band edge compared 
to the case of an inflection point. Namely, as follows from 1474f) in the case of a 
degenerate point of order 4 we have [fi, fi] = where in the case of an inflection point, 
according to Ij458jl . [fi,fi] 7^ 0. 

Using again (gHH> together with (j^ . (j4TT7|l and (jiT^I we get 

[fx, f 3 ] = [Tofi.Tofa] = [fi + f , f 3 + f 2 ] = [fi,fa] + [fo, fa] + [fi,f 2 ] + [fo, fa] , (475) 
readily implying 

[fo,fa] + [fi,f 2 ]=0. (476) 
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Summarizing P% )l -l|17T )l . l|17^ . P7I|) . g7<3 

[fo, fo] - 0, [f , fi] - 0, [fx.fi] = 0, [f , f 2 ] = 0, (477) 

and 

[fi,fa] + [fa,fi] =0, [fo,f 3 ] + [fi,f 2 ] = 0. (478) 
Notice that the first identity in (|478|l implies that [f i , f 2] is pure imaginary, i.e. 

Re{[fi,f 2 ]} = 0, [h,f 2 ] = ilm[h,f 2 ] . (479) 

In particular, the first three identities in l|477|l imply that 

for any f G Span{f , fj : [f, f] = 0. (480) 

As we have already pointed out this behavior of fluxes reflected by l|48U|) is very 
different from the case of an inflection point for which always [fi, fi] 7^ 0. 

9. 3. Transfer matrix and the fluxes for a degenerate point of order 2 

In the case of a degenerate point of order 2 the transfer matrix has a Jordan block of 
rank 2 and according to (|269[) there exists a basis fj, j = 0, 1, 2, 3 in C 4 for which we 
have 

Tof = Cofo, %h = Cofi + Cofo, ICol - 1, (481) 

%k = Cfa, C + Co, (482) 

where |C| = 1 or |£| 7^ 1. There are some additional relations not given here. Notice 
that for (|481|l the relation l|432[l applies yielding 

[fo,fo] = 0. (483) 

Observe that since C 7^ Co, m both cases we have |C| = 1 or \C\ 7^ 1, in view of l|131(l 
and H133(l . then 

[fo,f 2 ]=0. (484) 

Notice that C 7^ Co implies CoC 7^ 1- Using l|43T|) together with l |4*5Tt and lEfg!^ we 
obtain 

[fi,f 2 ] = [T f 1 ,T f 2 ] = CoC[fi + fo,f2] = CoC[fx,fa], (485) 
implying, in view of Co C 7^ 1, 

[fi,fa]=0. (486) 
Notice that (|438|l applied in this case yielding 

[fo,fi] + [fi,fo] = 0. (487) 

Evidently, 

Re [f , fi] = 0, [f , fi] = i Im [f , fx] . (488) 
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10. Perturbation theory for the matrix of reflection coefficients. 

In Section 6 we have introduce and studied the matrix p of reflection coefficients and 
its relation to the space St (0;w,k T ). In this section we study the behavior of the 
matrix p at frequencies u> close to the frequency of a degenerate point ujq, i.e. as 
v = uj — u — ► 0. To do that we first describe the space S T (0; u>, k T ) by the formula 
(|227(l where the vectors $i and $2 depend on the frequency v. As we know by now 
that this dependence has the form (see l|3U9|) and the section on the perturbation 
theory, and also l|317|l ) 

$j {v) = + + +0 (v 2 ) ,j = l, 2, (489) 
v = —\v~ = oiqv™ + O (v^J 

where n = 2,3,4 is the degeneracy index. To get an expansion for p(v) we use the 
relations ^525 ^ -^555 |l . First we need obtain an expansion for the matrices Q ± (v) 

Q ± (f) =Qv+Qfv + (v 2 ) , where Q± = Q ± (0) , (490) 
based on (|229|l and (|232j) . Notice that according to <)229|) we have 



3 = 1,2, 




) , «o = 


n! 


U (n) (fc )J 



Qo ± = Q ± (0) = [^ + ^ 2 + ] t [*i(0)$ 2 (0)] 1 



Z±$i(0)) (^,$2(0) 

^,$1(0)) ^±,$2(0) 



Pw,k T 

(491) 

We assume that the vectors $1 (z/) and $2 (f) are chosen so that for v = they are 
linearly independent, i.e. 

{$1 (0) , $2 (0)} = {$10, $20} are linearly independent. (492) 

The fulfillment of the condition l|492l) allows the limit space Span-j^ (v) , $2 as 
v — > to be described as the two-dimensional space Span{$i (0) , $2 (0)}- It is also 
necessary for the invertibility of the matrix Q ± (0) defined by l|491(l . i.e. for 

det Q ± (0) ^ 0. (493) 

In fact, we should always have 

det Q ± (u) ^ for any v, (494) 

for any semi-infinite slab problem. 

Based on the above, we get the following asymptotic expansions for [Q + (v)] 
and p (i>) 

[Q + M] ~* = [Qo] ~* - QoQt [Qo] ^t + O (v 2 ) , (495) 
p {u) = Q- (u) [Q+ (p)] _1 = {Q - + Q- X v + O {v 2 ) }{Qt + Qtv + (v 2 )}' 1 
= { Qv + Q^v + O (v 2 ) } { [Qt] - QtQt [Qt] ^t + O (v 2 ) } 

= Qo [QtV 1 - [Q0T 1 {QoQI -QoQt} [<#] _1 * + o(* a ) , 



Slow light in photonic crystals 



83 



or 



p (v) = p + p x v + O {v 2 ) , where (496) 

pa = Qo [QoV 1 > pi = [QoV 1 {QoQl - QtQt] [QtY 1 ■ 

The relation l|49t)[l readily implies 

P f 0) P (y) = pIpo + p\piv + pIpo" + O (/> 2 ) (497) 



r2 ( a+ . „) = k^j! + 2Re{( PlQ +, Po a+)} , + Q 2 
|a+| |a+| 

= r (a+;0) + — -g v + O (u ) 

\a+\ 

The relation (|498[> together with (|236|l yield the following expression for the flux 
associated with incident wave described by a + 

[$ (a + ;^) ,$ (a + ;u)] (499) 

= (1 - r 2 («+; 0) |« + | 2 = U - - - + O (^)) . 

\ |a + l Ia + I / 



11. Relevant modes near a degenerate point. 

In this section we describe in detail the properties of the space of relevant modes 
S t (0;lu) — <S T (0; a;, k r ) (suppressing in the notation its dependence on k T ) in a 
vicinity of a degenerate point for all the three cases, namely, an inflection point, 
n = 3, and band edges of orders n — 2,4. 

The general framework determining the basic properties of the space <S T (0; u) 
has been considered in the subsections "Basic properties of the space of relevant 
eigenmodes" and "Matrix of reflection coefficients and the flux quadratic form" . At 
this point having investigated the spectral properties of the transfer matrix T (v) , 
v = to — luq at a degenerate point luq and as v — > (see Section "Spectral perturbation 
theory of the transfer matrix a point of degeneracy"), we can provide more details 
of the properties of S T (0;cj) including the asymptotic behavior of the flux and the 
reflection coefficients of the relevant eigenmodes for a semi-infinite slab as v — > 0. 



11.1. Relevant modes near an inflection point. 

In this section we study the basic properties of the relevant eigenmodes and, in 
particular, the space S t (0;lo) as functions of the frequency lj in a vicinity of an 
inflection point ljq, i.e. for u> = luq + v when v is small. 

Using the equalities 1)3571) and (|450II - H454|I we get the following asymptotic 
formulae as v — > for the fluxes 



[eo (v) , eo (v)] - v 2 [fi,fi] - 2z> 2 Re [f , f 2 ] + O (v 3 ) = 

= 'ii? [fj.fi] + O (v 3 ) = 3a 2 oV 2 / 3 [fi, fi] + O (v) 



(500) 
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where, in view of Q339[) - (|34ip . 

. = a ,V3 )Hl = a 3 = _6_. (5Q1) 

Notice now that the relations Ijl311133|i . together with (|357H and (|426l) yield the 
following formulae for the fluxes 

[ei (v) , t x (i/)] = [ea (i/) , e 2 (i/)] = 0, [e 3 (i/) , e 3 (i/)] - [f 3l f 3 ] + O (i/) (502) 
[e (i>) , ei (v)] = [t (v) , e 2 [v)\ = [e (v) , e 3 (i/)] = 0, 
[e 3 iy) , e {v)\ = [e 3 (i/) , e x (i/)] = [e 3 (f) , e 2 (i/)] = 0. 

To handle in a uniform fashion both positive and negative v we introduce 

e» (V) = e ■ (v) - I e+ M if " " ° ^ M - I 9l W if " - (503) 
<iW-^W-( e _ (j) ) if V<Q > "i W _ \ fl 2 (i/) if ^<o ■ iaudj 

e+ (i/) = ei (i/) , e_ (i/) = e 2 (i/) . 
Then it follows from (JH32J>, dED that 

| = e -^l*l (l + o(>/ 3 )). (504) 

Notice that, as follows from (JSUUJ, the vector eo {v) has a positive flux and, hence, 
corresponds to a propagating mode. As to z\ [y), in view of l)504[l . it corresponds to 
an evanescent mode decaying as x 3 — > oo. So, based on H193[l . we obtain 

S T (0; uj + v) = Span {e (f) , e} (i/)} . (505) 
Then using l|361|) one verifies that the following limit exists 

S T (0; qjq) — lim S T (0; w + v) — Span {f , fi} . (506) 

v— >0 

Observe that the representation l|506|) for the space S T (0;wo) together with l|458|) and 
yield 

[Pofo + /3ifi, /3 fo + /?ifi] = l/3i| 2 [fi, fi] ^ if & ^ 0. (507) 

The relation (|507|l combined with l|22()[l imply the following very important property 
of the reflection coefficient r 0) at the inflection point uiq, i.e. v = 0, 

for almost all a + € C 2 : the reflection coefficient r (a + ; 0) < 1. (508) 

The relation l.5fl<S|) clearly indicates that the reflection coefficients r (a + ;0) are always 
strictly less than 1 /or aZZ i/ie relevant eigenmodes of the semi-infinite periodic stack, 
with the only exception when EM field value of the eigenmode at the surface of the slab 
is fo . In other words, at an inflection point there always will be a positive fraction of 
the incident energy transmitted through the infinite slab. In fact, by proper design of 
the slab one can achieve almost 100% transmission of the incident energy. In contrast, 
at any band edges the transmission is always exactly zero and the reflection is always 
100%, as we will see from the analysis in the following sections. 

More elaborate analysis yields asymptotic expressions for the matrix of reflection 
coefficients p, as determined by (|232[1 and l|496|) . and other related quantities for 
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nonzero by small v = to — ui^. Indeed, let us use in the relations i|229[l - i|232|) the 
vectors <&i and $2 defined by 



if. I \ l \ if. I \ ui I \ ( V ) ~~ e ° ( V ) 



(509) 



ll/ ^sign(iy) _ 1) 

Notice that (|356[) and (|359[) yield the following representation 

e (y) = (1 + ir 2 *> + O (i> 2 )) fo + (i* + ir^ 2 ) fi - ^ 2 f 2 + O (v 3 ) , (510) 

f)«(^) = [r 2 + 0(^)]fo + 



^ — sign(^) 1 ✓ 



i (?Mgu(6) - l) 

6 



1 l?sign(£) ~ J- J 



3C tl ' ~" " u w'" (fco) 
In particular, l|509l) and (|510|) yield for z/ = 

*i(Q)=fo, *i(0)=r 2 f +fi, 

implying that 



(511) 
(512) 



{$1 (0) , $2 (0)} are linearly independent. 

The relation l|512[l implies that the condition (|492[1 is satisfied. 

Now, let find the value $ (a + ; v) of the eigenmode corresponding to the incident 
wave a + . Using 1)228(1 . (|231[1 . (|495|l we consequently obtain 



¥>2 0) 



<Pi (0) 

f2 (0) 



0{y), 



<Pi (0) 

<P2 (0) 



[Q + (0)] '«+ (513) 



(514) 

= ^ (l/) $! (!/) + ^ M $2 M - Vl M e (I/) + ^ M % ( " } 1 ^ ! _ C " ' " ! 



(?sign(£) _ l) 



^2 (0) 



(^sign(w) 



-jy [ e Bign(^) (f) - eo (f)] + O (1) . 



Observe that the decomposition l|514|) of the vector <I> (a + ; v) into a linear combination 
of eigenvectors eo {v ) and e s ign(i>) (z/) of the transfer matrix T (y) signifies that 



the amplitude of the eigenmode inside the slab is — 



¥>2 (0) 



+ 0(1). (515) 



W (?sign(0) - l) 

Combining l|500|l , l|502|l with l|510|l , l|513|) we get the following formula for the flux 

[$ (a+;u) (a + ;u)] (516) 



?+ — 1 3 



The formula l|51t)[) together with (|237|l yield 

|[g+ (o)]- 1 ^} 2 

f (a+; 1/) = 1 - r 2 (a+; „) = i J2 [f l5 f x ] + O (t>) 

3 \a+\ 



(517) 
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More accurate computation based on (|499|l . ((496(1 and H5U8|I implies the following 
asymptotic formulae for the transmission and reflection coefficients 

,2 / + \ , 11 + \ , \Po a+ \ 2 2Re{(p 1 a+,p q+)} , , 2 . 

t (a + ;^) = f-r (a + ;^=f-—— 2 — 2 1/ + O (i/ ) , (518) 

|a+| |a+| 

2 ,.\ |p a + | 2 2Re{(p 1 a+,p a+)} _ / , 2 , |p a + | 2 „ 
r (<V, z/) = - — —g- H — —J f + C [y ) , 2 < 1, 



Po 



\ a I |« T | |OH , 

Qo [<#]~\ Pi= [QoT'l^Qr-Qo+Qi"} [QoT 1 - 



Observe that the formulae (|518|l involve the matrix p 1 requiring more terms in the 
expressions for the eigenvectors eo (v) and e± (v) (namely we need to compute the 
matrix A 2 as defined in the Appendix 2). When the exact value of the matrix p l is 

found we can find the exact value of the coefficient 2Rc {( p ^^ p ° a 11 i n ((518(1 . At this 

point we are interested in the concrete value of 2Rc i( Pl " ' p ° a )} an( j f or that reason 
we have not carried out the computation of the matrix A 2 . 

11.2. Relevant modes near a degeneracy point of order 4- 

In this section we study the basic properties of the relevant eigenmodes and, in 
particular, the space iSp (0; u) as functions of the frequency u> in a vicinity of a 
degenerate point loo of order n = 4, i.e. for to = ujq + v when v is small. Without loss 
of generality we assume v > 0. 

Notice that the eigenvector eo (v) corresponds to the eigenvalue 9\ (v) for which 

\6 Q {v)\ = l, (519) 

and, hence, the corresponding eigenmodc is a propagating one. 

Using the equalities 1(388(1 and 1(477(1 - 1(479(1 we get the following asymptotic 
formulae as v — > for the fluxes 

[e (y) , e (v)] = Im [fi , f 2 ] + O (z> 4 ) , (520) 
Pi] 3 i Pi] 42 , 1/4 4 4! 



8C ti' Z 8C ti' u ' ° wWfa)" 

Notice now in view of 1)384(1 the eigenvector ei (j/) corresponds to the eigenvalue d\ (y) 
for which evidently 

\e 1 {v)\ = |e ifeo+i ^| (l + 0(i) 2 )) = e~ 6 , v>0, (521) 

and, hence, the corresponding eigenmode is an evanescent one. 

In view of 1)519(1 and 1(521(1 . we can use the relations 1(131(1 and ((132(1 yielding 

[ ei (i/),ei (v)] = [e x (p) ,c (u)} = 0. (522) 

Hence, as follows from ((519(1 and 1(520(1 the vectors eo {v) and ei (z/) correspond 
respectively to a propagating and evanescent modes. So based on 1(193|) we get 

5 T (0; w + v) = Span {e (i/) , e x (i/)} . (523) 
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Then using i|393|) one verifies that the following limit exists 

5 T (0;w ) = limS T (0;w o + f) = Span {f , fi} . (524) 

v— >0 

Notice the representation (|524[1 for the space S T (0;wo) together with 148U|I yield 

[f,f] =0for anyf e«S x (0;o;o). (525) 

The relation H525|) combined with (|223[1 and (|224H implies the following very important 
property of the reflection coefficient r (a + ; 0) at the inflection point ojo, i.e. v = 0, 

for every a + € C 2 the reflection coefficient r = 1 and p = I 2 . (526) 

The relation $526)) clearly indicates that the reflection coefficient r(a + ;0) is always 
exactly 1 for all the relevant eigenmodes of semi-infinite periodic stack In other words, 
at any degenerate point of order n = 4, 100% of the incident energy is reflected and, 
hence, no energy is transmitted. In contrast, at any inflection point a positive fraction 
of the incident energy is always transmitted. 

More elaborate analysis provides asymptotic expressions for the matrix of 
reflection coefficients p, as determined by l|232|) and (|496|l . and other related quantities 
for nonzero, but small, v = u> — uiq. Indeed, let us use in the relations (|229[l - i|232|) the 
vectors $i and $2 defined by 

$! <y) = to (v) , $2 (u) = f)r („) = ^t^). . (527) 

Notice that (|356[) - H359[) yields the following representation 

e (v) - (1 + W 2 v + O {v 2 )) f + {iv + ir^ 2 ) fj - z/ 2 f 2 + O (v 3 ) , (528) 
Ijj, (u) = (rr + 3r 2 ) f + fi + O 0) . 

In particular, lf577|) and (|52~g|) yield for v = 

$1 (0) = f , $1 (0) = (ti + 3t 2 ) f + fi, 

implying that 

{$1 (0) , $2 (0)} are linearly independent. (529) 

The relation l|529[) shows that condition (|492ll is satisfied. 

Now let find the value $ (a + ; v) of the eigenmode corresponding to the incident 
wave a + . Using (|228H . (|231|l and (|495|l we consequently obtain 

$ (a + ; v) = 



•P2 W) 



<Pi (o) 

<P2 (0) 



0(y) 



<Pi (0) 

¥>2 (0) 



= [Q + (o)\ 



(530) 



= <Pi {v) $1 (v) + tp 2 (v) $2 (y) = <p t (v) e (y) + f 2 W) 



(531) 



gp (^) - ei (y) 
0(i+l) 



y 2 (Q) 



[eo(i/)-ei (i/)]+0(l). 



Slow light in photonic crystals 88 

Observe that the decomposition l|531|) of the vector $ (a + ; v) into a linear combination 
of eigenvectors eo (y) and t\ (is) of the transfer matrix T (y) signifies that 

the amplitude of the eigenmode inside the slab is f 2 ^\ +0(1). (532) 

Combining (|520|l and i|522[l with (|530|l and i|520[l we get the following formula for the 
flux 



i + 1 



2 
2 

The formula f53l)l together with (f2~3Tj) yields 
a /„.+ . ,.\ 1 „2 



V2 (0) 2 



2£Im[fi,f 2 ]+0(£ 2 ) (533) 



{[Q + (0)] 'a+^ImlhJ^ + O^ 2 



i 2 (a + ;^) = 1-r 2 (a + ;u) (534) 
{[Q+ (Or'a+I'lmlf!,^^ 



la+l 2 



-v + 0(v 2 ) 



An alternative computation based on (|499|l . (|496|l and (|526|l (which turns into 
! + | 2 = |a + | 2 for all a + ) impliei 
transmission and reflection coefficients 



p a + \ 2 = |<^ + | 2 for all ce + ) implies the following asymptotic formulae for the 



.11 + \ -, 2 / + \ 2Re{(p 1 a+,p a+)} , ,w,2\ ,„r^ 
t (a T ;^J=l — r (cr;i/) = ^ v + O [y ) , (535) 

lrv+ I 



-2\ 



r2 = t + 2^^+^+)} , + Q , 2 

|a+| 

p = Qo [Oo] _1 , Pi = [QoT* Wo Ql - Qo~Qi~} W]" 1 • 

J_Z.5. Relevant modes near a degenerate point of order 2 

In this section we study the basic properties of the relevant eigenmodes and, in 
particular, the space 5 T (0; u>) as functions of the frequency w in a vicinity of a 
degenerate point u>o of order n — 2, i.e. for u> = ujq + v when v is small. Without loss 
of generality we assume v > 0. 

Notice that the eigenvector eo (v) corresponds to the eigenvalue 9\ (v) for which 

6 (y) = e lfco+i * (l + O (v 2 )) : |0 O = 1, (536) 

with the corresponding eigenmode propagating in the positive direction. 

Using the equalities (|415|) . I|483|) and H486|) - (|488|) we get the following asymptotic 
formulae as v — > for the fluxes 



[eo (f) , eo (f)] = -2Im[fo,fi]i> + 0(^ 2 ) , £ = a a y/v, a = J (537) 

y (fc ) 

In particular, for ^ = the equality (|537|) implies 

[eo(0),eo(0)] = 0. (538) 
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Notice that in view of Q384[l the eigenvector ei (y) corresponds to the eigenvalue 9\ (v) 
for which 

Q {y) = e ik °-^ (l + O (v 2 )) : \6 X (y)\ = 1, v > 0, (539) 

with the corresponding eigenmode propagating in the negative direction. 

In view of l|536|) and 1539|) . we can use the relations l|131|) and (|133|) yielding 

[ei (v) , e (i/)] = 0. (540) 

So, unlike in situations for n — 3,4 in the case n = 2 only the vector eo (v) belongs 
to the space S T (0;u>o + v), where another one, namely t\{v), does not belong to 
<S T (0; u>o + v) since it corresponds to an eigenmode propagating in the negative 
direction. Hence, the second vector in S T (0; loq + v) must be one of e 2 (v) and e3 (v). 
Without loss of generality we assume that is e 2 (v), and, hence 

S T {0;u + is) = Span{e (i>),e 2 (i/)}. (541) 

Now there can be two possibilities: 1 6*2 {v)\ < 1 or |# 2 {v)\ — 1. The most interesting 
case is 

\6 2 { V )\ <1, v> 0, (542) 

when the corresponding mode is an evanescent one. Since we are interested to know 
if there can be any transmission of energy by a mode related to a regular band edge 
under assumption Ij542(l , the only possibility of the transmission will the mode related 
to the band edge. 

In the case of |#2 (z')! = 1 the corresponding mode will be a common one 
propagating in the positive direction with non-zero velocity. In this case the calculation 
is similar to the case (|542() with the only difference that we have to pick the single 
eigenmode related to the band edge and find the corresponding flux and the reflection 
coefficient. For that mode the result is the same as in the case (|542H . 

So, we assume now that the condition l|542|) is satisfied. Notice that under the 
condition (|542|) in view of l|131|) and i|132|) we have 

[e {y) , e 2 {v)\ = 0, [e 2 (u) , e 2 {v)\ - 0, v > (543) 

Then it follows from (|541|l that 

S T (0; too) = Span {e (0) , e 2 (0)} . (544) 

Notice the representation (|544|l for the space <S T (0;wo) together with l|538|) and l|543|l 
yields 

[f,f] =0for anyf e«S T (0;w ). (545) 

The relation l|545(l combined with (|223|l and (|224|l implies the following very important 
property of the reflection coefficient r (a + ; 0) at the inflection point loq, i.e. v = 0, 

for every a + € C 2 the reflection coefficient r (a + ) = 1 and p = 7 2 . (546) 

The relation \52b)) clearly indicates that the reflection coefficient r (a + ; 0) is always 
exactly 1 for all the relevant eigenmodes of semi-infinite periodic stack related to the 
band edge. 

More elaborate analysis yields asymptotic expressions for the matrix of reflection 
coefficients p, as determined by (|232|l and l|496|) . and other related quantities for 
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nonzero, but small, v = u> — uj a . Indeed, let us use in the relations i|229|) - (|232[l the 
vectors $ x and <& 2 defined by 

$1 (v) = e (v) , $2 0) = e 2 (i/) . (547) 

Notice that generically eo (0) and e 2 (0) are always linearly independent and, hence, 

{$i (0) , $2 (0)} are linearly independent. (548) 

The relation l|548[) shows that condition (|492ll is satisfied. 

Now let us find the value <!> (a + ; v) of the eigenmode corresponding to the incident 
wave a + . Using (|228[) . ()23ip . (|495|l we consequently obtain 

$ (a + ;u) = 



<P2 M 



<Pi (o) 

<P2 (0) 



0{v), 



<Pi (o) 

¥>2 (0) 



[Q + (o) 



(549) 



$ (a+; i/) = ip x (y) $i (i/) + <^ 2 (f) $2 (v) = <p x (0) e (0) + ip 2 (0) e 2 (0) + O iy) (550) 

Observe that the decomposition l|550|) of the vector <I> (a + ; 1/) into a linear combination 
of eigenvectors eo (y) and ei (y) of the transfer matrix T [y] signifies that 

the amplitude of the eigenmode inside the slab is ip 1 (0) + O (y) . (551) 

Combining l|537|) , l|543[) with l|549|) , l|550|) we get the following formula for the flux 

[$ (a+; 1/) , $ (a+; 1/)] = -2 | Vl (0)| 2 Im [f , fi] !> + O (^ 2 ) (552) 

= -2 { [Q+ (0)] _1 Im [fo, fx] * + O (v 2 ) . 

The formula together with (|257j) yields 

2{[Q+ (Or'a+flm [fo, fi] 

t 2 («+; 1/) = 1 - r 2 (a+; = ^ + O (*> 2 ) . (553) 

|a+| 

An alternative computation based on (|499|l . (|496|l and (|546|l (which turns into 
|/5 Q! + | = |a + | for all a + ) implies the following asymptotic formulae for the 
transmission and reflection coefficients 

.2/ + \ 1 2 / + \ 2Re{(p 1 q + ,p a+)} , „ / , 2 \ / KK/n 

£ \pr;v) = l-r z {oT\v) = v + 0[y), (554) 

\a+\ 

2t + \ , , 2Re{(p 1 a+,p a+)} , , /2 , 

r (a + ; 1/) = 1 + — - 2 v + O [v ) , 

\a+\ 

Po = Qo [QoT\ pi = [QoV 1 {QoQI - QoQt} [Q0T 1 ■ 

11.4- Asymptotic analysis summary 

The final results on the reflection coefficients, transmission and flux are formulated in 
Section "relevant modes at degenerate points" . 
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We reiterate that for any relevant eigenmode of a semi-infinite slab the following 
fundamental relation holds for its energy flux [<f>, $] and the reflection and transmission 
coefficients 



t 2 (a+) = 1 - 



[$ (a+), $(a+)] 



a" 1 



where the two-dimensional vector a + describes the incident wave in a properly chosen 
basis and $ (a + ) is the corresponding EM field at the surface of the slab. 

One the most important quantitative results of the analysis of the reflection 
coefficient in a vicinity of band edges and inflection points is summarized by the 
following formulae for the reflection coefficient r &s v = uj — cjq^O 



inflection point n = 3 : r 2 



sign v 



regular band edge n = 2 : r 2 = 1 — Cq\u 



degenerate band edge n = 4 : r 2 = 1 — cq\v 



| 1/3 , 0<r o <l; 
1/2 / v > for a lower edge 
' \ v < for an upper edge 

1/4 J v > for a lower edge 
i/ < for an upper edge 



where c s i gnu denotes one of the constants c± corresponding to the sign of v. 

The above formulae for the reflection coefficient indicate clearly that on approach 
to a band edge the reflection coefficient approaches 1. In contrast, in a vicinity of an 
inflection point the reflection coefficient approaches a less than 1 number r and can 
be made arbitrarily small for properly designed structures. 

The table below shows the asymptotic behavior of the slow mode group velocity, 
the saturation amplitude, and the semi-infinite slab transmittance as the frequency 
approaches the respective stationary point. 



Rank of degeneracy 


Group velocity 


Saturation amplitude 


Transmittance 


2 (regular band edge) 


1 1 1/2 


1 




3 (stationary inflection pint) 




1 1 — 1/3 

M 




4 (degenerate band edge) 


wr 


1 1-1/4 

M 





Using this table we summarize the basic properties of the eigenmodes at 
frequencies close to the band edges and inflection points as follows. 

1. For a regular band edge there are no energy relevant Floquet modes, and as 
v — > the group velocity and the maximal flux vanish as M 1 ^ 2 , whereas the saturation 
amplitude remains finite. The light does slow down in the vicinity of a regular band 
edge, but only a vanishing fraction enters the photonic slab, while the rest is reflected 
back to space. 

2. For a stationary inflection point there is a relevant Floquet mode at loq, and 
as v = uj — loq — > the group velocity vanishes at the rate |^| 2 ^ 3 with the maximal 
flux remaining finite, and the saturation amplitude diverging as |^| -1 ^ 3 . The slab 
transmittance at u> = loq remains finite and can even be close to unity. 

3. For a 4-fold degenerate band edge, there is a relevant non-Bloch Floquet 
mode. As v — > 0, the respective slow mode group velocity vanishes as |^| 3 ^ 4 , while the 
saturation amplitude diverges as |^p 1//4 . The transmitted energy flux, along with the 
slab transmittance, vanishes as \v^^ A ' . 
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12. Summary 

Although the existence of slow electromagnetic modes in photonic crystals is quite 
obvious, the next question is whether and how such modes can be excited by incident 
light. In other words, we need to know whether or not a significant fraction of the 
incident light energy can be converted into a slow mode with virtually zero group 
velocity in a semi-infinite photonic crystal. We have shown that it can be done, 
but only if the slow mode is associated with a stationary inflection point of the 
dispersion relation oj{k). In this special case, the incident light with the proper 
frequency, polarization, and direction of incidence is completely converted into the 
slow frozen mode with huge amplitude and vanishingly small group velocity. Such a 
phenomenon constitute the frozen mode regime. By contrast, if a slow electromagnetic 
mode relates to a photonic band edge, the incident wave will be reflected back to space. 
Not every photonic crystal can have the dispersion relation with a stationary inflection 
point and, thereby, support the frozen mode regime. For instance, one-dimensional 
periodic arrays (periodic layered structures) must include specially oriented anisotropic 
layers, in order to support the proper dispersion relation and the frozen mode 
regime. Photonic crystals with three-dimensional periodicity are not required to have 
anisotropic constitutive components. 

Generally, the possibility of conversion of an incident wave into a slow mode 
appears to be directly related to the character of the respective Bloch dispersion 
relation ui (k) of the periodic structure. This fundamental relation exists regardless of 
the specific physical realization of the periodic structure supporting such a dispersion 
relation. For instance, although neither periodically modulated waveguides, nor 
periodic arrays of coupled resonators are formally photonic crystals, still, as soon as the 
respective Bloch electromagnetic dispersion relation develops a stationary inflection 
point, one can expect the frozen mode regime at the respective frequency. 

Not every periodic array can have the dispersion relation with a stationary 
inflection point. Symmetry-based considerations similar to those applied above to 
the case of periodic layered arrays, can provide a meaningful guidance on how to find 
the proper structure. 
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Appendix 1: basic properties of the transfer matrix 

According to (|91|l . the 4x4 matrix T (v), v — uj — ojq satisfies the following 
identity 

T- 1 {v) = JT* (u) J. (555) 
The identity implies, in particular that 

|detT(i/)|=l (556) 

if C is an eigenvalues of T (v) then l/£* is also an eigenvalue. (557) 

In other words, the statement (|557|l yields that if ( = pe 1 ^ is the polar form of an 
eigenvalue Q and p ^ 1 then 

C = pe ict > and C = p^e 10 are both eigenvalues of T (v) (558) 
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and that if £ is an eigenvalue of T (y) then l/£* is an eigenvalue too. 

The above properties of eigenvalues of T (y) imply the following statements. 

1. Suppose that (,{v) is an eigenvalue of T (y) depending on v continuously. 
Suppose also that £ (0) has multiplicity 1 and |C(0)| = 1. Then there exists a 
sufficiently small <5 > such that 



|C(»| = 1 for any \v\ < 5. 
To show Q559JI we need the following elementary implication: 



if |C | = land K-Col <e< 2 then 



C* 



Co 



< 6e. 



Assume now that the assumption of (|560[) holds. Notice that 

i _. c 



and 



icr 



||C|-l| = HC|-|Coll<lC-Col<e<-. 



with then implies that 

Now using (|561(l we get 
1 

C 
< 



1 — e < |C| < 1 + e. 



Co 



icr 



-Cc 



icr 



-1 C + C-Co 



(559) 
(560) 

(561) 

(562) 
(563) 

(564) 



+ IC-Col 



1 - icr 

ICI 



+ ic-Coi 



The inequalities l|562|l . (|564|l together with e < \ yield 



C* 



Co 



l-(l-e) 



3-2e n 
1 — e 



(565) 



which is the desired inequality (|560ll . 

Using the fact the C (0) has multiplicity one and applying the standard 
perturbation theory arguments we can always find < eo < 1 an d Sq > such 
that 

for \v\ < Sq the eigenvalue C {v) is the only one in the circle |£ — £ (0)| < eo- (566) 
Now using the continuity of £ (v) we can always find a positive S < Sq such that 

(567) 

u u 

Observe that (|565|l and (|567|l imply that 



for \v\<5< S we have |C M - C (0)1 < — < -. 

6 6 



C(o) 



< e for \v\< 5 < S . 



(568) 
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Assume for the sake of the argument that for some \v\ < 5 we have |£ (v)\ ^ 1. Then 
based on l(568|) and general properties of T (v) we have to conclude that ^wn/, 7^ C { v ) 
is another eigenvalue of T (v) residing in the circle |£ — ( (0)| < eo. But this clearly 
contradicts to ((566(1 implying the desired relation (|559|l . 

2. Suppose that for < \v\ < S the matrix T (y) has four different eigenvalues 
QAv), j = 1, ...4 each continuously depending on v and having the following 
properties: 

There exists a Co sucn that lim C (y) = Co, 3 = 1, 2, 3, and C4 (0) ^ Co- (569) 

In other words, for small \v\ the eigenvalue C4 (^) has multiplicity one and is well 
separated from the other 3 different eigenvalues Cj (v), j — 1,2,3 which converge as 
v —* to a Co- 
Then we claim that 

|Col =1 and |C 4 (0)1 = 1, (570) 
and there exists a S > such that 

Kj = ^ ^ or a ^ l eas t one J = 1, 2, 3 and |f| < <5. (571) 



|C 4 WI = lfor |i/|<<5. 



(572) 



Observe, first, that the first equality in (|569|l follows from the conditions of 1)569(1 and 
general properties of T (y), since if | Co I 7^ 1, we would have three more eigenvalues 
,7 = 1,2,3 with the total number of eigenvalues 6. That is, of course, impossible 



for 4x4 matrix implying that the first equality in (|569|l holds. As to the second, it 
follows from proving that |Col = 1 an d the identity Q556f) . since then we must have 

ICol 3 IC 4 (0)1 = 1. 



To show l|571f) we use the limit conditions in (|569|l and the fact that 



CM 



is also 



an eigenvalue. Indeed, if, for the sake of the argument, we assume that (|569|l does not 
hold we have to conclude that in a infinitesimally small vicinity of Co there will be at 
least four different eigenvalues which is impossible in view of the second condition in 
l|569|) . This completes the proof of 1)571(1 . As to the proof of ((572|) it follows from the 
statement 1(559(1 . 

Appendix 2: perturbation theory for a diagonal matrix 

Particular constructions of the perturbation theory we discuss here follow from 
HO] and g7j. Suppose that 



W (v) = W + vWi + v 2 W 2 + 
and Tq is a diagonal matrix with distinct elements, i.e. 

Wl ... 
w 2 ' ■ : 



(573) 



Wo 







'■■ 

W n 



where w m 7^ Wj if m^J. 



(574) 
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To diagonalize T (v) we use the approach outlined in |46| and used in |47| . Namely, 
there exists the following representation for the diagonal form X of T (v) 

e -S{u) We S{u) =Wq + v ^ + ^ +... } S(v)=vS! + v 2 S 2 + ■■■ , 



C 

where the matrices Si, S 2 , . . . do not depend on v and Xj_, X 2 , 
find Sj and Xj we use the Hausdorf 's representation 



, (575) 
are diagonal. To 



X = e~ s We s = W + [W, S] 



[[W,S\,S\ 



where 



(576) 



where the brackets denote the commutator of two matrices 

[A, B]=AB- BA. 

Substituting (|575|) into (|576J) and equating the terms of like powers in v, we obtain 
the following expressions for the matrices Xj 

X x = [T , Si] + Ti , X 2 = [T , S 2 ] + T 2 + [T x , S x ) + ~ [[T , S x ] , S x ] , . . . . (577) 

To find Xj we introduce for a matrix Y its representation as the sum of its diagonal 
Y part and the remaining part Y with zero diagonal elements (so called integrable 
matrix |4fi| ) 

Y = Y + Y, y = diag(T), y = y-diag(F). (578) 

Then 

X 1 = \W ,Si]+Wi + W 1 , (579) 

and to get rid of the integrable part W\ of W\ we take Si to be the solution of the 
equation 

[W ,S 1 ] = -W 1 . 

The solution to this equation is 

1 



[Si 



Wi 



J I" 



3 ± m; [Si] .. = 0. 



Consequently, 

X-l =diag(W r 1 ) = Wi. 
To find X2 we recast the equation (157711 as 

X 2 = [W ,S 2 ] +Y 2 + Y 2 , Y 2 = W 2 + [WuSi] + X - [[W 0l Si] , Si] 
Applying to this equation the same approach as for l|579|) we get 



X 2 = diag (Y 2 ) = Y 2 , Y 2 = W 2 + [Wi,Si] + - [[Wo, Si] , Si] 



(580) 
(581) 
(582) 

(583) 
(584) 



[S-2 



jrn 



Y 2 



jrn 



j ^ m; [S 2 ]^ = 



Using HH3-ISH3 



recast l|584|) 



X 2 = diag W 2 



[S 2 



1 

z 2 



Wi,Si 
, j 7^ to; [S 2 ] 



(585) 



0; Z 2 



W 2 



Wi - -Wi,Si 



Slow light in photonic crystals 



96 



[1] L. Brillouin. Wave Propagation and Group Velocity. (Academic, New York, 1960). 
[2] L. D. Landau, E. M. Lifshitz, L. P. Pitaevskii. Electrodynamics of continuous media. (Pergamon, 
N.Y. 1984). 

[3] A. Yariv and Pochi Yeh. Optical Waves in Crystals. ("A Wiley-Interscience publication", 1984). 
[4] A. Sommerfeld, Phys. Z. 8, 841 (1907). 

[5] A. Kuzmich, A. Dogariu, L. J. Wang, P. W. Milonni, R. Y. Chiao, Phys. Rev. Lett. 86, 3925 

(2001) . 

[6] R. W. Boyd, D. J. Gauthicr, in Progress in Optics, E. Wolf, Ed. (Elsevier, Amsterdam, 2002), 
vol. 43. 

[7] P. W. Milonni, J. Phys. B 35, R31 (2002). 

[8] V. G. Veselago. The electrodynamics of substances with simultaneously negative values of e and 

IX Soviet Physics USPEKHI 10, 509-514 (1968). 
[9] L. Hau, S. Harris, Z. Dutton, C. Bchroozi. Light speed reduction to 17 metres per second in an 

ultracold atomic gas. Nature, 397, 594 - 598, (1999). 
[10] M. Kash, V. Sautcnkov, Al. Zibrov, L. Hollberg, G. Welch, M. Lukin, Yu. Rostovtsev, E. Fry, and 

M. Scully. Ultraslow Group Velocity and Enhanced Nonlinear Optical Effects in a Coherently 

Driven Hot Atomic Gas. Phys. Rev. Lett. 82, #26, 5229 - 5232 (1999) 
[11] D. Budker, D. F. Kimball, S. M. Rochester, and V. V. Yashchuk. Nonlinear Magneto-optics and 

Reduced Group Velocity of Light in Atomic Vapor with Slow Ground State Relaxation. Phys. 

Rev. Lett. 83, # 9, 1767 (1999). 
[12] M. Lukin and A. Imamoglu. Controlling photons using clcctromagnctically induced transparency. 

Nature, 413, 273 - 276, (2001). 
[13] D. F. Phillips, A. Fleischhauer, A. Mair, and R. L. Walsworth, M. D. Lukin. Storage of Light 

in Atomic Vapor. Phys. Rev. Lett. 86, # 5, 783 - 786, (2001). 
[14] A.V. Turukhin, V. S. Sudarshanam, M. S. Shahriar, J. A. Musscr, B. S. Ham, and P. R. Hemmer. 

Observation of Ultraslow and Stored Light Pulses in a Solid. Phys. Rev. Lett. 88, # 2, 023602 

(2002) . 

[15] Matthew S. Bigelow, Nick N. Lcpcshkin, and Robert W. Boyd. Observation of Ultraslow Light 
Propagation in a Ruby Crystal at Room Temperature. Phys. Rev. Lett. 90, # 11, 113903 

(2003) . 

[16] Matthew S. Bigelow,* Nick N. Lepeshkin, Robert W. Boyd. Superluminal and Slow Light 

Propagation in a Room-Temperature Solid. Science, 301, 200 (2003). 
[17] M. Yanik and S. Fan. Stopping Light All Optically. Phys. Rev. Lett. 92, # 8, 083901, (2004). 
[18] J. Heebner and R. Boyd. Slow and stopped light. 'Slow' and 'fast' light in resonator-coupled 

waveguides. Journal of modern optics, 49, #14/15, 2629 (2002). 
[19] J. Heebner and R. Boyd. Slow light, induced dispersion, enhanced nonlinearity and optical 

solitons in a resonator-array waveguide. Phys. Rev. E65, 036619 (2002). 
[20] A. Mclloni, F. Morichetti, M. Maritclli. Linear and nonlinear pulse propagation in coupled 

resonator slow-wave optical structures. Optical and Quantum Electronics 35, 365 (2003) 
[21] J. Poon, J. Scheuer, Y. Xu, and A. Yariv. Designing coupled- resonator optical waveguide delay 

lines. J. Opt. Soc. Am. B, Vol. 21, No. 9 (2004). 
[22] J. Scheuer, G. Paloczi, J. Poon, and A. Yariv. Toward the Slowing and Storage of Light. OPN, 

16, 36 (2005). 

[23] J. B. Khurgin. Optical buffers based on slow light in clcctromagnctically induced transparent 
media and coupled resonator structures: comparative analysis. J. Opt. Soc. Am. B 22, #5, 
1062(2005). 

[24] J. B. Khurgin. Expanding the bandwidth of slow-light photonic devices based on coupled 

resonators. Optic Letters, 30, # 5, 513 (2005). 
[25] M. Notomi,l K. Yamada,2 A. Shinya,l J. Takahashi,2 C. Takahashi,2 and I. Yokohama. 

ExtremeJy Large Group-Velocity Dispersion of Line-Defect Waveguides in Photonic Crystal 

Slabs. Phys. Rev. Lett. 87, #25, 253902 (2001). 
[26] M. Scalora, R. J. Flynn, S. B. Reinhardt, R. L. Fork, M. J. Bloemer, M. D. Tocci, C. M. Bowden, 

H. S. Ledbetter, J. M. Bendickson, J. P. Dowling, R. P. Leavitt. Ultrashort pulse propagation 

at the photonic band edge: Large tunable group delay with minimal distortion and loss. Phys. 

Rev. E54, #2, R1078 (1996). 
[27] M. Bloemer, K. Myneni, M. Centini, M. Scalora, and G. D'Aguanno. Transit time of optical 

pulses propagating through a finite length medium. Phys. Rev. E65, 056615 (2002). 
[28] M. Soljacic, S. Johnson, S. Fan, M. Ibanescu, E. Ippen, and J. D. Joannopoulos. Photonic- 

crystal slow-light enhancement of nonlinear phase sensitivity. J. Opt. Soc. Am. B., 19, #9, 

2052 (2002). 

[29] A. Figotin, and I. Vitebskiy. Electromagnetic unidircctionality in magnetic photonic crystals. 



Slow light in photonic crystals 



97 



Phys. Rev. B67, 165210 (2003). 
[30] A. Figotin, and I. Vitebskiy. Oblique frozen modes in layered media. Phys. Rev. E68, 036609 
(2003). 

[31] J. Ballato, A. Ballato, A. Figotin, and I. Vitebskiy. Frozen light in periodic stacks of anisotropic 

layers. Phys. Rev. E71, (2005). 
[32] S. Molchanov and B. Vainberg. Slowdown of the wave packages in finite slabs of periodic media. 

Waves Random Media 14, 411 (2004) 
[33] A. Figotin and I. Vitebskiy. Gigantic transmission band-edge resonance in periodic stacks of 

anisotropic layers. Phys. Rev. E72, 036619, (2005). 
[34] A. Figotin and I. Vitebskiy. Electromagnetic unidircctionality and frozen modes in magnetic 

photonic crystals. JMMM, 300, 117 (2006). 
[35] S. E. Harris. Electromagnctically induced transparency. Physics Today 50, # 7, 36 (1997). 
[36] J. Joannopoulos, R. Meade, and J. Winn. Photonic Crystals. (Princeton University Press, 1995). 
[37] Pochi Ych. "Optical Waves in Layered Media", (Wiley, New York, 1988). 

[38] Wcng Cho Chew. " Waves and Fields in Inhomogencous Media" , (Van Nostrand Rcinhold, New 
York, 1990). 

[39] M. Notomi. Theory of light propagation in strongly modulated photonic crystals: Rcfractionlikc 
behavior in the vicinity of the photonic band gap. Phys. Rev. B62, 10696 (2000) 

[40] A. Figotin, and I. Vitebsky. Nonreciprocal magnetic photonic crystals. Phys. Rev. E63, 066609 
(2001). 

[41] D. W. Berreman. J. Opt. Soc. Am. A62, 502-10 (1972). 

[42] I. Abdulhalim. Analytic propagation matrix method for anisotropic magneto-optic layered 

media, J. Opt. A: Pure Appl. Opt. 2, 557 (2000). 
[43] I. Abdulhalim. Analytic propagation matrix method for linear optics of arbitrary biaxial layered 

media, J. Opt. A: Pure Appl. Opt. 1, 646 (1999). 
[44] R. Bellman. Introduction to Matrix Analysis. (SIAM. Philadelphia, 1997) 

[45] E. Coddington and R. Carlson. Linear Ordinary Differential Equations. (SIAM, Philadelphia, 
1997). 

[46] V. N. Bogacvski and A. Povzner, Algebraic Methods in Nonlinear Perturbation Theory, Springer- 
Verlag, New York, 1991. 

[47] Figotin A. and Godin Yu. Spetral Properties of Thin-Film Photonic Crystals, SIAM J. APPL. 

MATH., 2001, Vol. 61, No. 6, pp. 1959-1979. 
[48] Kato T., Perturbation Theory of Linear Operators, Springer, 1995. 
[49] Lancaster P. and Tismcnctsky M., The Theory of Matrices, Academic Press, 1985. 
[50] Wilkinson, J., The Algebraic Eigenvalue Problem, Oxford University Press, 1996. 



